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Abstract 

This  report  describes  the  use  of  computer  program  ’SEEK’  which  works  in 
conjunction  with  two  user-written  subroutines  and  an  input  data  file  to 
perform  an  optimization  procedure  on  a  user’s  problem.  The  optimization 
method  uses  a  modified  feasible  directions  gradient  technique.  SEEK  is 
written  in  ANSI  standard  Fortran  77,  has  an  object  size  of  about  46  K  bytes 
and  can  be  used  on  a  personal  computer  running  DOS.  This  report  describes  the 
use  of  the  program  and  discusses  the  optimizing  method.  The  program  use  is 
illustrated  with  four  example  problems:  a  bushing  design,  a  helical  coil 
spring  design,  a  gear  mesh  design  and  a  two-parameter  Wei  bull  life-reliability 
curve  fit. 
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SUMMARY 


This  report  describes  the  use  of  a  computer  program,  ’SEEK’  for 
engineering  design  optimization.  The  program  is  not  complete  in  itself  in 
that  it  is  written  to  work  with  problem  specific  user  subroutines  and  an  input 
data  file.  It  performs  a  gradient  search  optimization  of  the  user’s  problem 
to  find  an  optimal  set  of  design  parameter  values.  Optimization  is  performed 
using  a  modified  feasible  directions  gradient  technique.  The  program  is 
written  in  ANSI  standard  Fortran  77  and  has  an  object  size  of  about  46  K  bytes 
for  the  optimizing  code  alone  for  use  on  a  personal  computer  running  DOS.  Its 
source  code  is  about  1,200  lines  in  length  and  its  size  is  39  K  bytes. 

In  the  OPTIMIZATION  FORMAT  section,  the  four  interface  vectors  to  the 
procedure  are  described.  These  vectors  are:  1)  the  problem  constants,  2)  the 
independent  design  parameters,  3)  the  constraint  bounds,  and  4)  the  objective 
function  terms.  The  problem  constants  and  independent  design  parameter  values 
define  a  specific  trial  design.  In  the  optimization,  the  program  varies  the 
design  parameter  values  to  search  for  the  design  which  has  the  best  objective 
function  value  while  satisfying  the  constraint  bounds. 

In  the  PROGRAMMING  section,  the  report  describes  the  two  analysis 
routines,  BOUNDS  and  VALUES  which  the  user  must  write  to  evaluate  the 
constrained  functions  and  the  objective  function.  BOUNDS  takes  as  input  the 
constant  and  design  parameter  vectors  and  produces  a  vector  of  the  constrained 
function  values  as  its  output.  VALUES  takes  as  input  the  same  constant  and 
design  parameter  vectors  and  produces  the  objective  function  value  vector  as 
its  output.  The  format  and  make-up  of  the  input  data  file  which  gives  the 
four  vectors  of  constants,  design  parameters,  constraint  limits  and  objective 
function  weighting  coefficients  are  described.  In  this  file,  extensive 
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labeling  for  the  four  vectors  is  included.  The  design  verification  feature  of 
the  optimization  program  is  described  also.  Once  a  numerical  optimum  is 
found,  the  program  provides  the  user  the  opportunity  to  try  alternate  designs 

for  comparison  purposes. 

To  illustrate  the  use  of  the  program,  four  examples  are  presented: 

1)  a  bushing  design,  2)  a  helical  coil  spring  design,  3)  a  spur  gear  mesh 
design,  and  4)  a  Weibull  data  curve  fit.  The  bushing  design  problem  is  to 
find  the  length  and  diameter  of  a  bushing  to  minimize  the  friction  torque  in 
the  bearing  for  a  given  load  and  material  properties.  The  spring  design  is  to 
find  the  wire  diameter  and  mean  coil  diameter  which  support  a  given 
alternating  load  with  a  required  stiffness.  Three  different  design  objectives 
of  minimum  spring  weight,  minimum  spring  height  and  minimum  coil  volume  are 
sought.  The  gear  design  is  to  determine  the  number  of  pinion  teeth,  diametral 
pitch  and  face  width  for  a  compact  spur  gear  set  to  transmit  a  given  power  at 
a  given  input  speed  with  a  given  speed  reduction.  Two  objectives  are  sought: 
1)  minimum  center  distance  for  a  desired  life  and  2)  maximum  life  for  a  given 
center  distance.  The  fourth  example  shows  the  fitting  of  a  two-parameter 

Weibull  distribution  to  life  data  for  a  series  of  identical  units  tested  at 

the  same  load. 

Each  example  includes  an  analysis  of  the  problem,  the  organization  of 

the  optimization  variables,  the  writing  of  the  analysis  subroutines  and  the 

input  data  file  for  a  specific  problem,  the  compiling  and  running  of  the 
program,  the  use  of  design  verification  to  obtain  reasonable  values  and  an 
interpretation  of  the  output  data  file.  In  each  case,  the  design  verification 
opportunity  leads  to  the  discovery  of  a  practical  solution  with  near  optimal 
characteristics. 


2 


After  the  examples,  the  OPTIMIZATION  METHOD  section  presents  a 
description  of  the  gradient  search  method  with  its  three  modes  of  operation: 

I)  unconstrained  searching,  2)  acceptable  design  searching,  and  3)  feasible 
direction  searching  along  one  or  more  design  constraints.  The  unconstrained 
searching  mode  employs  the  gradient  in  the  objective  function  to  improve  the 
design  at  the  fastest  rate  when  possible.  The  acceptable  design  searching 
uses  the  sum  of  the  gradients  in  the  violated  constraints  to  find  the 
acceptable  design  region.  And,  the  feasible  direction  searching  mode  combines 
the  gradients  in  the  objective  function  and  the  violated  constraints  to 
improve  the  objective  function  value  while  keeping  the  trials  within  the 
acceptable  design  region.  This  section  concludes  with  a  description  of  the 
program  structure  and  operation. 
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INTRODUCTION 


Optimization  is  a  mathematical  process  of  seeking  the  most  favorable 
combination  of  parameters  to  achieve  the  best  outcome  possible  [1-3].  In 
design,  one  constantly  searches  for  an  ideal  trade-off  of  conflicting 
performance  objectives.  For  aircraft  transmissions,  for  example,  we  might 
wish  to  obtain  the  lightest  transmission  which  has  some  minimum  acceptable 
service  life,  or  we  may  want  to  maximize  the  service  life  at  a  given 
transmission  weight  [4].  These  objectives  are  sought  throughout  the  design 
and  development  process  with  repetitive  design  descriptions  and  evaluations 
on  paper  or  CAD  layout  at  the  design  stage  and  in  hardware  at  the  prototype 
stage.  Many  of  the  optimizing,  trade-off  decisions  which  develop  and  improve 
the  product  are  made  by  engineers  without  the  help  of  mathematical  models  of 
the  product’s  performance. 

Optimization  offers  the  promise  of  assistance  with  the  difficult  trade¬ 
off  decisions  at  the  early  stages  of  the  design,  before  costly  prototypes  are 
constructed  and  tested  [5].  A  spectrum  of  optimization  codes  have  been 
written  to  assist  in  the  design  of  complex  structures  which  can  be  modeled 
with  large  matrices  of  simultaneous  linear  equations  [2,3].  With  an  objective 
computer  search  through  the  space  of  potential  designs,  we  are  given  a  greater 

chance  of  determining  a  truly  optimum  design. 

Unfortunately,  many  mechanical  design  problems  cannot  be  described  by  a 
set  of  linear  equations.  The  size  of  the  problem  as  measured  by  the  number  of 
independent  design  parameters  and  the  constraints  which  are  to  be  satisfied 
may  be  small.  But  the  complexity  of  the  problem  may  be  large  due  to  the  non¬ 
linear  and  often  discontinuous  nature  of  the  design  property  and  constraint 
relations.  This  combination  keeps  the  optimization  of  mechanical  designs  in 
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the  realm  of  an  art  which  requires  considerable  engineering  insight  to 
complement  the  available  mathematical  models  and  computer  solutions  [6]. 

In  modeling  mechanical  systems,  one  is  often  confronted  with  the  choice 
of  obtaining  exact  solutions  to  an  approximate  model  or  obtaining  approximate 
solutions  to  an  exact  model.  With  the  designer  active  in  the  process,  rapid 
solutions  of  an  approximate  model  can  lead  to  practical  optimal  designs  even 
when  the  mathematical  optimum  contains  some  flaws  [7].  Modified  gradient 
optimization  techniques  such  as  the  feasible  directions  method  are  quite 
powerful  and  rapid  when  they  have  continuous  models  in  which  gradients  can  be 
calculated.  The  method  of  this  optimization  is  to  have  the  engineer  write 
Fortran  subroutines  which  model  the  design  with  continuous  properties  as 
functions  of  the  constants  and  independent  design  parameters  which  define  the 
design.  The  optimizer  can  then  find  the  optimal  solution  to  this  ’ideal’ 
problem  and  report  it  to  the  engineer  to  allow  a  check  of  alternate  designs 
which  satisfy  additional  constraints  of  practicality.  The  end  result  is  a 
practical,  optimal  design. 

This  report  describes  the  use  and  background  of  a  Fortran  program,  SEEK, 
which  is  written  to  assist  the  mechanical  designer  in  developing  balanced 
optimal  designs.  The  intent  of  the  program  is  to  keep  the  engineer  in  the 
process  while  providing  a  systematic  search  of  potential  designs.  In  doing 
this,  it  allows  the  engineer  to  use  the  mathematical  models  of  the 
optimization  to  evaluate  near  optimal,  practical  designs.  SEEK,  which 
requires  two  user-written  analysis  subroutines,  reads  its  input  from  an  ASCII 
data  file  and  writes  the  output  of  the  optimization  both  to  the  screen  and  to 
an  output  ASCII  log  file.  To  document  the  optimization  clearly,  the  input 
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data  file  includes  a  significant  amount  of  text  to  describe  the  numerical 
values  in  the  output  file. 

The  report  includes  an  OPTIMIZATION  FORMAT  section  which  describes  the 
basic  format  of  the  optimization  problem:  the  constants,  design  parameters, 
constraint  bounds  and  objective  function;  as  well  as  a  PROGRAMMING  section, 
which  describes  the  programming  required  to  prepare  SEEK  for  use  in  an 
interactive  design  session.  This  section  describes:  the  two  analysis 
subroutines,  the  input  data  file,  the  use  of  the  program  with  changes  to  the 
input  file  and  design  verification  in  the  interactive  session. 

To  demonstrate  the  power  and  ease  of  use  of  this  optimization  procedure, 
several  small  design  examples  follow  in  the  next  sections:  a  bushing,  a 
spring,  a  spur  gear  mesh,  and  the  fitting  of  a  two  parameter  Weibull 

distribution  to  experimental  test  data. 

An  OPTIMIZATION  METHOD  section  follows  which  describes  the  structure  and 
operation  of  the  optimization  code.  This  code  is  small  with  1200  lines  and 
less  than  40  k  bytes  so  the  optimization  can  be  performed  on  a  personal 
computer  running  with  DOS.  The  speed  of  a  486  machine  may  become  attractive 
for  the  more  complex  analysis  models. 
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OPTIMIZATION  FORMAT 


An  optimization  problem  may  be  formulated  as  a  constrained  search  in 
terms  of  four  vectors  and  two  sets  of  relations.  In  this  formulation,  only 
inequality  relations  are  used  for  the  constraints.  The  four  vectors  are: 

1)  the  constants  of  the  problem  which  do  not  change  for  a  given  design,  2)  the 
parameters  which  define  a  design  and  which  are  the  sought  values,  3)  the 
constraint  values  which  may  be  upper  or  lower  bounds  on  properties  of  the 
design,  and  4)  the  objective  function’s  weighting  coefficients. 

In  this  formulation,  at  least  two  Fortran  subroutines  are  needed: 

1)  BOUNDS  which  evaluates  the  constrained  variable  values  in  terms  of  the 
constants  and  design  parameter  values  and  2)  VALUES  which  evaluates  the 
objective  function’s  value  for  a  given  set  of  constants  and  design  parameter 
values.  These  two  subroutines  combine  with  the  input  data  file  to  define  the 
specific  problem  for  optimization.  The  gradient  calculations  which  perform 
the  optimization  by  calling  these  two  subroutines  repeatedly  and  the  input  and 
output  file  interfacing  are  contained  in  the  SEEK  Fortran  code. 

Constants 

Each  problem  is  defined  by  a  series  of  constant  values  such  as:  size, 
power  level,  speed  of  operation,  elastic  modulus,  material  strength  or 
requested  service  life.  These  constants  are  fixed  for  all  trial  designs  of 
the  optimization  effort,  and  the  constrained  properties  and  objective  function 
values  are  direct  functions  of  them.  The  constants  may  change  for  a  different 
design  using  the  same  analysis  subroutines,  however.  For  example,  designs 
made  of  steel  will  have  different  properties  from  those  made  of  nylon,  yet 
each  steel  design  will  have  the  same  material  stiffness  and  strength  as  the 
other  steel  designs. 
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The  program  will  read  these  constants  and  their  labels  from  the  ASCII 
input  file,  store  them  in  arrays  and  use  the  values  whenever  the  constrained 
property  or  objective  function  values  are  calculated. 

Parameters 

In  each  problem,  we  are  searching  for  a  set  of  parameter  values  which 
optimize  the  objective  function  to  either  a  minimum  or  a  maximum  value.  These 
parameters  are  the  second  vector  entered  in  the  input  data  file.  The 
optimization  scheme  proceeds  by  analyzing  repeated  trials  until  it  selects  one 
for  which  the  analysis  yields  an  optimum  objective  function  value.  So  the 
values  entered  for  the  design  parameters  include  an  initial  value  for  each 
parameter  for  the  first  trial  analysis.  This  vector  also  includes  upper  and 
lower  values  for  each  design  parameter.  These  upper  and  lower  values  serve  to 
establish  the  relative  sensitivity  of  the  parameter  for  the  gradient  searches, 
but  do  not  limit  the  value  itself.  By  increasing  the  span  between  the  upper 
and  lower  values  for  a  design  parameter,  the  user  can  increase  the  sensitivity 
of  that  parameter  in  the  design  search.  If  it  appears  that  a  design  parameter 
is  not  changing  as  the  optimizer  seeks  out  better  designs,  increasing  this 
span  between  the  upper  and  lower  values  will  increase  its  tendency  to  change 
in  future  optimization  runs. 

After  reading  these  parameter  values  and  their  labels  from  the  input 
file,  the  program  will  store  them  in  arrays,  use  the  parameter  ranges  to  set 
relative  sensitivities  which  will  not  change  throughout  the  optimization  and 
place  the  initial  parameter  values  into  the  parameter  array.  The  parameter 
array  will  change  throughout  the  operation  of  the  program  until  it  contains 
the  values  of  the  parameters  which  optimize  the  objective  function. 
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Constraints 


Limiting  each  design  is  a  series  of  constraints  on  the  properties  of  the 
design.  These  constraints  may  be  applied  directly  to  one  or  more  of  the 
parameters  such  as  the  thickness  of  a  beam  or  they  may  be  applied  to  a 
calculated  property  of  the  design  such  as  the  maximum  stress  in  the  beam.  The 
constraints  of  this  algorithm  are  inequality  constraints.  In  the  general 
optimization  formulation,  two  types  of  constraints  are  possible:  inequality 
and  equality.  Gradient  search  algorithms  require  a  continuum  of  parameter  and 
property  values  in  which  to  move  around  in  search  of  the  optimum.  Inequality 
constraints  provide  boundaries  to  the  design  space  but  do  not  diminish  it. 

Equality  constraints  reduce  the  space  by  one  dimension.  Each  equality 
constraint  transforms  one  independent  design  parameter  into  a  dependent  design 
parameter.  There  are  two  ways  to  include  an  equality  constraint  in  this 
algorithm:  1)  reduce  the  design  space  by  one,  or  2)  enter  the  equality 
constraint  into  the  data  as  an  inequality  constraint. 

The  first  method  is  the  best  because  it  simplifies  the  calculations 
making  the  optimization  more  direct  and  faster.  If  the  width  of  a  rectangular 
spring  is  always  seven  times  its  height,  then  one  can  remove  the  width  from 
the  list  of  independent  design  parameters  and  set  it  equal  to  seven  times  the 
height  in  the  calculations.  This  reduces  all  gradient  calculations  by  one 
element  and  leaves  a  full  design  space  for  the  remaining  parameters’  gradient 
searches. 

The  second  method  may  be  easier  to  implement  if  the  equality  constraint 
is  not  tied  to  the  parameters  directly.  By  entering  it  as  an  inequality 
constraint,  one  leaves  the  design  space  at  its  larger  dimensional  size  and 
cuts  it  in  half  with  the  bounding  value.  Since  the  unbounded  optimal  design 
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would  probably  lie  off  this  constraint,  placing  this  bound  between  the 
unbounded  optimum  and  the  acceptable  design  space  will  place  the  bounded 
optimal  design  right  on  this  constraint  and  thus  actually  satisfy  the  equality 
constraint. 

The  program  will  read  these  constraints,  their  directions  and  their 
labels  from  the  input  file,  place  them  in  arrays  and  use  them  in  the 
subroutine  BOUNDS,  which  is  provided  by  the  user,  to  limit  the  acceptable 
design  space  throughout  the  optimization. 

Objective  Function 

In  each  optimization,  some  property  or  combination  of  properties,  called 
the  objective  function,  is  to  be  minimized  or  maximized.  The  weighting 
coefficients  of  these  properties  are  the  last  vector  entered  in  the  input  data 
file.  These  coefficients  may  be  percentages,  unit  conversions  to  place  the 
properties  in  the  same  dimensions  or  they  may  be  switches  such  as  zero  and  one 
to  change  the  optimization  in  the  data  file  by  changing  the  sought  objective. 
The  assumption  is  that  the  objective  function  to  be  optimized  may  be  expressed 
as  a  linear  sum  of  terms,  each  with  its  own  weighting  coefficient.  The 
weighting  coefficients  and  direction  of  optimization  will  remain  fixed 
throughout  the  optimization. 

The  program  will  read  these  coefficients  and  their  labels  from  the  end 
of  the  input  data  file,  place  them  in  arrays  and  use  the  coefficients  to 
modify  the  objective  function  property  values.  These  values  are  calculated  in 
subroutine  VALUES  which  is  provided  by  the  user. 
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PROGRAMMING 


With  these  four  vectors  defined  and  labeled,  the  program  starts  the 
output  file  with  an  echo  of  the  input  data  to  document  the  optimization 
problem.  It  then  proceeds  to  seek  an  optimal  design  with  the  modified 
feasible  directions  gradient  method.  Using  gradients  in  the  objective 
function  and  in  the  violated  constraints,  the  program  can  move  from  an  initial 
design  which  does  not  satisfy  the  design  constraints  to  designs  which  are 
valid.  It  can  also  improve  a  valid  initial  design  to  obtain  an  optimum  design 
within  the  assumptions  of  the  model. 

Once  an  optimum  design  is  found,  the  program  prints:  the  found  design 
parameters,  the  objective  function  value  with  its  component  function  values, 
and  the  constraints  with  both  their  design  and  limit  values.  The  program  then 
offers  the  user  the  opportunity  to  try  alternate  designs.  On  receiving  the 
revised  design  parameter  values,  the  program  uses  subroutines  BOUNDS  and 
VALUES  to  check  this  design,  prints  out  its  properties  and  offers  the  user  the 
chance  to  try  another  alternate  design.  All  design  trials  are  printed  to  the 
screen  and  the  ASCII  output  log  file. 

Analysis  Subroutines 

To  model  a  problem,  the  program  needs  two  analysis  subroutines:  BOUNDS 
and  VALUES.  These  subroutines  are  problem  specific  and  should  match  the  input 
data.  Subroutine  BOUNDS  calculates  the  constrained  property  values  for  each 
design  trial  and  gradient  perturbation  as  direct  functions  of  the  constants 
and  design  parameters  only.  Data  are  passed  to  BOUNDS  with  two  dynamically 
dimensioned  arrays:  CONST(NCO)  for  the  constants  and  X(NX)  for  the  design 
parameters,  and  the  constraint  value  results  are  returned  to  the  program  in 
the  array  VCSTR(NCS).  Subroutine  VALUES  calculates  the  objective  function 
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values  also  as  direct  functions  of  the  constants  and  design  parameters  only. 
Data  are  passed  to  VALUES  with  the  same  two  arrays:  CONST(NCO)  for  the 
constants  and  X(NX)  for  the  design  parameters,  and  the  object  function  values 

are  returned  in  the  array  OBJECT(NOB). 

Both  subroutines  must  be  able  to  determine  their  outputs  as  continuous 
functions  over  the  range  of  design  parameters  used.  Since  small  perturbations 
are  given  to  the  design  parameters  to  determine  corresponding  changes  in  the 
objective  function  value  and  in  the  constrained  variable  values  throughout  the 
design  search,  the  subroutine  calculations  must  be  defined  and  continuous. 
Discrete  parameter  requirements  such  as  integer  tooth  numbers  for  gears  and 
standard  component  sizes  can  be  added  by  the  user  in  the  verification  stage  of 
the  optimization  process.  They  cannot  be  included  in  the  simulation  model 
itsel f . 

These  subroutines  may  contain  formulas,  interpolated  data,  iterations  or 
other  subroutines  as  long  as  the  resulting  calculations  yield  continuous 
functions  of  the  design  parameters.  If  the  subroutines  call  other 
subroutines,  they  should  not  have  the  same  names  as  those  subroutines  included 
in  the  optimizing  part  of  SEEK,  which  are  listed  in  Table  1. 

Common  blocks  may  be  used  by  the  subroutines,  but  SEEK  uses  four  common 
blocks,  which  should  not  be  altered:  CURVE,  PAR,  VAR,  and  UNITS.  One  of  these 
common  blocks,  UNITS,  contains  four  integer  variables,  NW,  NR,  NF  and  ND. 

These  are  the  file  numbers  for  writing  and  reading  to  the  interface  devices. 

NW  identifies  the  screen,  NR  identifies  the  keyboard,  NF  identifies  the  output 
file  and  ND  identifies  the  input  file.  The  user  may  add  additional 
information  to  the  output  with  these  file  numbers,  bearing  in  mind  that  BOUNDS 
and  VALUES  are  called  many  times  by  SEEK  in  the  optimization  process. 
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Table  1 


Subroutines  of  Program  SEEK 


Line 

Name 

Function 

756 

BACK 

Search  for  Acceptable  Design  Space 

983 

BOUNCE 

Find  Gradient  Sum  of  the  Violated  Constraints 

1207 

BOUNDS 

User  Supplied  Constraint  Analysis 

1113 

CHECK 

Test  for  Constraint  Violation 

1033 

GRADNT 

Evaluate  a  Gradient 

1174 

MERIT 

Evaluate  the  Objective  Function  Sum 

724 

RESIZE 

Unscale  the  Design  Parameters  into  Real  Units 

912 

SCAN 

Increment  the  Design  in  the  Acceptable  Design 

Space 

874 

SCOUT 

Try  a  New  Design  Position  And  Check  the 

Constraints 

740 

SIZE 

Scale  the  Design  Parameters  to  Unit  Space 

1088 

UNIT 

Normalize  a  Vector 

1210+ 

VALUES 

User  Supplied  Objective  Function  Analysis 

1146 

WALL 

Evaluate  a  Specified  Constraint 
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Input  Data  File 

Coordinated  with  these  two  required  subroutines  is  the  ASCII  input  data 
file,  which  is  described  in  Appendix  A  -  SEEK.DOC.  The  initial  line  in  the 
data  file  is  a  text  line  of  fifty  characters  or  less  which  describes  the 
design  being  optimized.  This  is  followed  by  a  line  which  contains  a  single 
number,  NCO,  -  the  number  of  constant  values  to  follow,  which  is  the  first 
vector  in  the  data  file.  Each  optimization  constant  is  then  entered  in  a  set 
of  three  lines:  1)  the  numerical  value,  2)  the  name  of  the  constant  in  thirty 
characters  or  less,  and  3)  the  units  for  the  constant  in  twelve  characters  or 
less.  With  this  information,  the  program  will  label  the  constant  values 
whenever  it  prints  them. 

Following  the  constants  in  the  input  data  file  is  the  list  of 
independent  design  parameters,  which  is  the  second  vector.  After  the  last 
constant  has  been  listed,  the  next  line  is  once  again  a  single  number,  NX,  - 
the  number  of  parameter  values  to  follow.  Each  parameter  is  then  entered  in  a 
set  of  three  lines:  1)  three  numerical  values  -  a  low  estimate  for  the 
parameter,  a  high  estimate  and  an  initial  estimate;  2)  the  name  of  the 
parameter  in  thirty  characters  or  less;  and  3)  the  units  for  the  parameter  in 
twelve  characters  or  less. 

The  list  of  constraint  bounds  is  the  third  vector.  After  the  last 
parameter  has  been  listed,  the  next  line  is  a  single  number,  NCS,  -  the  number 
of  constraint  bounds  to  follow.  Each  bound  is  then  entered  in  a  set  of  three 
lines:  1)  the  word  ’UPPER’  or  ’LOWER’  followed  by  the  numerical  value 
including  its  decimal  point,  2) -the  name  of  the  constraint  in  thirty 
characters  or  less,  and  3)  the  units  of  the  constraint  in  twelve  characters  or 
less. 
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Finally,  the  list  of  weighting  coefficients  is  the  fourth  vector 
entered.  After  the  last  constraint  has  been  listed,  the  next  line  is  a  three 
letter  prefix,  ’MIN’  or  ’MAX’,  which  describes  the  direction  of  optimization. 
This  is  followed  by  a  line  with  a  single  number,  NOB,  -  the  number  of 
weighting  coefficients  to  follow.  Each  coefficient  is  then  entered  in  a  set 
of  three  lines:  1)  the  numerical  value,  2)  the  name  of  the  property  in  thirty 
characters  or  less,  and  3)  the  units  for  the  property  in  twelve  characters  or 
less. 

Design  Verification 

As  described  at  the  start  of  this  section,  the  use  of  SEEK  is  somewhat 
interactive.  Because  the  user  must  add  at  least  two  subroutines  to  the 
program  in  addition  to  the  input  data  file,  the  combined  program  must  be 
compiled  separately  for  each  optimization  application.  After  writing  the 
analysis  subroutines  with  the  proper  dynamic  dimensioning  in  the  calling 
arguments,  the  user  may  add  subroutines  to  the  end  of  the  source  code  for  SEEK 
and  compile  the  program  in  the  environment  in  which  it  is  to  be  run.  The 
compiler  should  be  a  Fortran  77  compiler  of  which  there  are  several  PC 
versions  available.  Once  compiled  and  linked  into  an  executable  program,  the 
optimizer  can  be  run  with  the  matching  data  file  to  find  an  optimal  design. 

With  the  interaction  of  the  data  file  and  the  analysis  subroutines,  the 
user  may  change  the  way  an  optimization  is  conducted  through  small  changes  in 
the  data  file.  By  using  one  and  zero  as  weighting  coefficients  to  an 
objective  function  that  contains  totally  different  terms  and  by  switching 
constraints  from  UPPER  to  LOWER  or  changing  their  values  to  make  them  active 
or  inactive,  one  can  change  an  optimization  significantly. 
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For  example,  one  could  have  a  transmission  life  optimization  program 
which  included  bounds  on  the  transmission  size  and  life  as  well  as  terms  in 
the  objective  function  for  size  and  life  [8].  By  requesting  that  the  size  be 
less  than  some  value,  that  the  life  be  greater  than  zero  and  by  having  a 
weighting  coefficient  of  one  for  life  and  zero  for  size  and  by  selecting  ’MAX’ 
in  the  input  data  file,  one  would  have  an  optimization  that  would  maximize  the 
life  of  the  transmission  within  a  given  acceptable  size.  Shifting  the 
requests  in  the  input  file  to  request  that  the  size  be  greater  than  zero,  that 
the  life  be  greater  than  some  desired  value,  and  by  having  a  weighting 
coefficient  of  zero  for  life  and  one  for  size  with  ’MIN’  selected  in  the  input 
data  file  would  minimize  the  size  of  the  transmission  for  the  requested 
service  life. 

Smaller  changes  in  limit  values  or  problem  constants  could  change  the 
size  of  a  requested  design  or  some  other  feature  without  requiring  a  change  in 
the  compiled  program.  As  stated  earlier,  the  program  generates  a  complete  log 
file  of  the  obtained  designs  and  the  verified  designs  in  response  to  keyboard 
input  after  an  ’ideal’  design  has  been  found  and  written  to  the  screen  and  the 
log  file.  Speed  of  execution  of  the  program  is  entirely  dependent  upon  the 
complexity  of  the  analysis  models.  Small  optimization  programs  can  run 
quickly  on  the  personal  computer.  The  following  four  sections  will 
demonstrate  the  use  of  SEEK  for  several  different  design  optimizations. 
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BUSHING  OPTIMIZATION 


Four  examples  of  increasing  complexity  will  be  presented  to  illustrate 
the  capability  of  SEEK.  The  first  example  is  that  of  the  design  of  a  low- 
speed  bearing  to  support  a  radial  shaft  load.  For  this  application,  the 
simplest  bearing  is  a  bushing  which  is  defined  by  its  material,  length  and 
diameter. 

Theory 

Consider  the  design  of  a  bushing  to  support  a  radial  shaft  load.  With 
little  or  no  lubrication,  a  bushing’s  capacity  is  both  strength  and  power 
limited  [9].  By  constraining  the  bearing  length  to  be  less  than  or  equal  to 
some  percentage  of  the  shaft  diameter,  one  can  treat  the  radial  load  as 
supported  uniformly  over  the  length  of  the  bushing.  Thus: 

L 

—  <  P  (1) 

D 


The  nominal  contact  pressure  in  the  bushing  can  then  be  taken  as: 


F 

P  =  -  (2) 

L  D 

where  the  pressure,  P,  is  measured  in  MPa;  the  load,  F,  is  in  Newtons;  and  the 
length,  L,  and  diameter,  D,  are  in  mm.  And  the  sliding  velocity  in  the 
bushing,  1^^,  measured  in  m/s  is: 


V 


s 


2  77 
60 


(3) 


where  the  shaft  speed,  Q,  is  in  RPM.  The  strength  limit  on  the  bushing  is 
thus: 
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The  power  limit  on  the  bushing,  which  is  proportional  by  the  coefficient  of 
friction  to  the  power  lost  in  the  bushing,  is  the  PV  factor  of: 


P  V, 


F 

L  D 


(5) 


P  V 


s 


F  Q 
L 


PV 


max 


(6) 


Figure  lisa  graph  of  the  contact  pressure  in  a  bushing  versus  the 
contact  sliding  velocity  which  shows  the  regions  of  the  two  pressure  limits. 
Acceptable  designs  have  pressures  lower  than  the  plotted  values.  As  the  speed 
increases,  the  power  limit  becomes  active  and  restricts  the  design  to  lower 
and  lower  acceptable  pressures.  Values  for  these  limits  for  both  metallic  and 
non-metallic  materials  are  readily  available  [10,11]. 

Given  adequate  strength,  a  bushing  may  be  sized  to  minimize  the 
frictional  torque  on  the  shaft.  This  torque  is  given  in  N-m  by: 


T 


f 


(7) 


Programming 

The  problem  of  designing  a  bushing  is  now  defined  mathematically.  The 
constants  which  specify  the  particular  application  are:  1)  the  radial  load,  F; 
2)  the  shaft  speed,  Q;  and  3)  the  coefficient  of  friction,  //.  The  two  design 
parameters  to  be  selected  are:  1)  the  bushing  length,  L;  and  2)  the  shaft 
diameter,  D.  The  three  inequality  constraints  are:  1)  the  length  to  diameter 
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ratio,  2)  the  acceptable  pressure,  and  3)  the  acceptable  pressure 

times  velocity  factor,  PV„,^.  All  three  constraints  are  upper  bounds.  The 
objective  function,  which  is  to  be  minimized,  is  the  frictional  torque,  T^. 
These  quantities  are  summarized  in  Table  2. 

The  relations  for  the  constraints  are  equations  (1),  (4)  and  (6),  and 
the  relation  for  the  objective  function  is  equation  (7).  A  subroutine  BOUNDS 
which  is  written  to  determine  the  constrained  values  using  equations  (1),  (2) 
and  (6)  is  listed  in  Table  3.  And  a  subroutine  VALUES  which  is  written  to 
determine  the  objective  function  value  using  equation  (7)  is  listed  in 
Table  4. 

The  simplicity  of  the  subroutines  matches  the  simplicity  of  the 
relations.  In  each  subroutine,  the  input  constant  and  design  parameter 
vectors  are  converted  to  individual  variables  which  have  names  that  identify 
them  more  clearly.  Then  the  equations  are  entered  in  an  easily  checked  format 
and  the  results  are  transferred  to  the  output  constrained  variable  or 
objective  function  vectors.  Note  that  the  vector  quantity  subscripts  match 
the  input  data  order.  These  quantities  are  numbered  in  the  output  file  echo 
of  the  input  to  assist  the  user  in  verifying  that  the  proper  input  and  output 
quantities  are  used  for  the  equations  in  the  two  analysis  subroutines. 

Once  written,  these  two  subroutines  must  be  compiled  and  linked  to 
program  SEEK  to  generate  an  executable  program  to  perform  the  optimization. 

One  way  to  do  this  is  to  add  the  subroutines  to  the  end  of  the  source  file  for 
program  SEEK. FOR,  save  the  combined  program  with  a  problem  specific  name  such 
as  BUSHING. FOR  and  compile  it.  The  result  will  be  an  executable  file, 
BUSHING.EXE.  A  second  way  would  be  to  compile  SEEK. FOR  and  the  two  analysis 
subroutines  BOUNDS. FOR  and  VALUES. FOR  separately  to  generate  object  files 
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Table  2 

Bushing  Optimization  Parameters 


Constants 

Design  Parameters 

Inequality  Constraints 

Objective  Function 

F 

L 

L 

-  <  B 

P  ^max 

^Vmin 

Q 

D 

P  <  P  ^ 

fJ 

max 

?  V  <  PV 
s  max 

f 
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Table  3 


SUBROUT I NE  BOUNDS ( CONST , NCO , X , NX , VCSTR , NCS ) 

BOUNDS  DETERMINES  THE  PRESENT  CONSTRAINT 
FUNCTION  VALUES 

FOR  A  BUSHING  DESIGN  EXAMPLE 

PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NCS  -  NUMBER  OF  INEQUALITY  CONSTRAINTS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  PARAMETERS 

VCSTR  -  PRESENT  CONSTRAINT  VALUES 

X  -  PRESENT  DESIGN  PARAMETER  VALUES 

ALL  VALUES  ARE  IN  PROBLEM  UNITS 

CONST(l)  =  F  -  RADIAL  LOAD  (POUNDS) 

CONST(2)  =  N  -  SHAFT  SPEED  (RPM) 

CONST(3)  =  f  -  FRICTION  COEFFICIENT 

X(l)  =  L  -  BUSHING  LENGTH  (IN) 

X(2)  =  D  -  BUSHING  DIAMETER  (IN) 

VCSTR(l)  =  P  -  AVERAGE  BUSHING  CONTACT  PRESSURE 

(PSI) 

VCSTR(2)  =  PV  -  BUSHING  PRESSURE  TIMES  VELOCITY 

FACTOR  (PSI  -  FT/MIN) 

VCSTR(3)  =  L/D  -  BUSHING  LENGTH  TO  DIAMETER  RATIO 


DIMENSION  CONST(NCO) ,X(NX) ,VCSTR(NCS) 
PI  =  3.14159265 
FORCE  =  CONST(l) 

RPM  =  CONST(2) 

BLEN  =  X(l) 

DIA  =  X(2) 

PRESS  =  FORCE/(BLEN*DIA) 

PV  =  0.001*PI*FORCE*RPM/(60.0*BLEN) 

RATIO  =  BLEN/DIA 

VCSTR(l)  =  PRESS 

VCSTR(2)  =  PV 

VCSTR(3)  =  RATIO 

RETURN 

END 
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Table  4 


Bushing  Objective  Function  Evaluation  Subroutine  Values 


SUBROUTINE  VALUES (CONST , NCO, X , NX , OBJECT, NOB) 

VALUES  DETERMINES  THE  PRESENT  DESIGN 
OBJECTIVE  FUNCTION  VALUES 

FOR  A  BUSHING  DESIGN  EXAMPLE 

PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NOB  -  NUMBER  OF  OBJECTIVE  FUNCTION  TERMS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  VARIABLES 

OBJECT  -  PRESENT  OBJECTIVE  FUNCTION  VALUES 

X  -  PRESENT  DESIGN  PARAMETER  VALUES 


ALL  VALUES  ARE  IN  PROBLEM  UNITS 


CONST(l)  =  F  -  RADIAL  LOAD  (POUNDS) 

CONST(2)  =  N  -  SHAFT  SPEED  (RPM) 

CONST(3)  =  f  -  FRICTION  COEFFICIENT 


X(l)  =  L  -  BUSHING  LENGTH  (IN) 

X(2)  =  D  -  BUSHING  DIAMETER  (IN) 


OBJECT(l)  =  Tf  -  BUSHING  FRICTION  TORQUE  (LB  -  IN)  C 


DIMENSION  CONST(NCO) ,X(NX) ,OBJECT(NOB) 
FORCE  =  CONST(l) 

FRICT  =  CONST(3) 

DIA  =  X(2) 

TORQUE  =  0.001*FRICT*F0RCE*DIA/2.0 

OBJECT(l)  =  TORQUE 

RETURN 

END 
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only.  These  object  files  can  then  be  linked  with  SEEK. OBJ  listed  first  to 
produce  an  executable  file,  BUSHING.EXE  [12,13]. 

Since  this  is  a  two  parameter  design  problem  with  a  single  objective 
function,  one  can  draw  two  graphs  which  illustrate  the  optimization.  The 
first  is  called  a  design  space  in  that  it  is  a  graph  in  coordinates  which 
match  the  design  parameters.  Points  in  the  graph  represent  specific  design 
parameter  values  or  designs.  Plotted  in  this  graph  are  the  design  constraint 
limits.  These  constraint  limits  divide  the  design  space  of  potential  designs 
into  two  regions:  1)  an  acceptable  design  region  in  which  all  design 
constraints  are  satisfied,  and  2)  an  unacceptable  design  region  in  which  at 
least  one  design  constraint  is  violated.  Figure  2  is  a  graph  of  a  design 
space  for  the  bushing  design  problem. 

The  second  graph,  Figure  3,  is  a  plot  of  the  objective  function  versus  a 
design  parameter.  If  the  objective  function  were  a  function  of  both  design 
parameters,  a  contour  plot  on  the  same  coordinates  as  Figure  2  would  be 
reguired  to  show  how  the  objective  function  varies  for  different  designs. 

Since  the  objective  function  of  equation  (7)  is  not  a  function  of  the  bushing 
length,  a  simple  graph  of  friction  torque  versus  shaft  diameter  shows  how  this 
property  varies  for  the  potential  designs.  Figure  3  is  drawn  directly  below 
the  design  space  so  that  the  objective  function  of  friction  torque  can  be 
visualized  as  a  plane  contour  rising  in  the  design  space  above. 

On  inspection  of  these  graphs,  it  is  obvious  that  the  optimal  design  is 
the  bushing  with  the  smallest  shaft  diameter  which  satisfies  the  three 
inequality  constraints  plotted  in  Figure  2.  Once  this  information  is  known, 
there  is  no  need  to  go  through  a  formal  optimization  to  find  the  optimal 
design.  Optimization  techniques  have  their  greatest  value  for  problems  for 
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which  optimum  solutions  are  not  yet  known.  Knowing  this  optimum  will  help  us 
verify  the  effectiveness  of  the  modified  gradient  method.  But  once  an  optimum 
solution  is  known,  either  from  a  graphical  analysis  or  by  a  computer 
optimization,  it  is  more  efficient  to  calculate  it  directly  [14]. 

Numerical  Example 

For  an  example,  consider  the  design  of  a  bushing  to  support  750  N  at  a 
shaft  speed  of  40  RPM.  The  shaft  is  steel  and  the  bushing  is  to  be  nylon 
which  has  a  coefficient  of  friction  with  steel  of  0.2  and  which  has  a  design 
pressure  limit  of  14  MPa  and  a  design  PV  limit  of  0.11  MPa  m/s.  In  this 
design,  the  length  is  to  be  limited  to  be  less  than  or  equal  to  seventy 
percent  of  the  shaft  diameter  and  the  shaft  size  and  bushing  length  are  to  be 
in  whole  mm’s. 

Table  5  is  a  listing  of  an  input  file  for  this  problem,  including  line 
numbers  which  are  not  part  of  the  input  file.  The  first  line  is  the  problem 
title.  The  next  line  is  the  number  of  constants,  3,  which  is  followed  by 
three  sets  of  three  lines.  Each  constant  is  identified  by  its  value,  its  name 
and  its  dimension. 

The  frictional  coefficient  has  no  dimension,  so  its  dimension  line  is 
left  blank  in  order  for  the  following  line  which  is  the  number  of  independent 
design  parameters,  2,  to  appear  in  its  proper  place.  If  this  line  11  were  not 
left  blank,  an  error  checking  routine  would  identify  an  error  in  the  input 
file  on  line  13  and  stop  the  program.  Line  13  is  the  next  line  of  text  in  the 
input  file  -  the  name  of  the  first  design  variable.  Due  to  the  missing 
line  11,  the  program  would  try  to  read  this  text  as  the  numerical  value  for 
the  first  design  variable.  Although  this  error  message  does  not  point 
directly  to  the  cause  of  the  reading  error,  it  does  indicate  the  presence  of 
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Table  5 


Line 

1 


First  Bushing  Input  File 


Column 

1 

RADIAL  NYLON  BUSHING 
3 

750.0 

RADIAL  LOAD 

NEWTONS 

40.0 

SHAFT  SPEED 
RPM 
0.2 

FRICTION  COEFFICIENT 
2 

0.0  10.0  5.0 
BUSHING  LENGTH 
mm 

0.0  10.0  5.0 

BUSHING  DIAMETER 
mm 
3 

UPPER  14.0 
CONTACT  PRESSURE 
MPa 

UPPER  0.11 
PV  FACTOR 
MPa  -  m/s 
UPPER  0.7 
LENGTH  TO  DIAMETER 
RATIO 
MIN 
1 

1.0 

FRICTION  TORQUE 
N  -  m 
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an  error. 


Checking  the  data  on  line  13  and  the  lines  that  precede  it  should 
lead  to  this  discovery  in  a  short  amount  of  time. 

The  next  six  lines  contain  the  design  parameter  values,  names  and 
dimensions.  The  value  lines  contain  three  numbers:  1)  the  low  estimate, 

2)  the  high  estimate,  and  3)  the  initial  estimate.  At  this  point  in  the 
solution,  we  know  the  least  about  the  values  to  enter  for  these  design 
parameters.  Let  us  guess  ranges  from  zero  to  ten  mm  and  initial  values  of 
five  mm  for  the  two  design  variables.  The  next  line  contains  the  number  of 
design  constraints,  3.  Following  are  nine  lines  with  the  three  constraint 
limit  types  and  values  with  their  decimal  points  on  the  first  lines,  their 
names  on  the  second  lines  and  their  units  on  the  third  lines. 

The  data  for  the  objective  function  vector  follows.  The  next  line 
contains  the  letters  ’MIN’  to  identify  minimization  as  the  direction  of 
optimization.  This  is  followed  by  a  line  with  the  single  value  of  one  to 
indicate  that  the  objective  function  has  only  one  term.  The  last  three  lines 
are  the  weighting  coefficient  value  and  the  name  and  units  for  the  term. 

This  file  is  saved  with  a  name  such  as  NYLONl.IN.  The  compiled  program 
BUSHING.EXE  can  now  be  run  by  typing  BUSHING  at  the  prompt.  As  shown  in 
Figure  4,  the  program  will  request  the  prefix  for  an  input  file,  which  should 
be  NYLONl  in  this  case.  Since  "data  points"  is  not  in  the  first  constant 
description,  the  program  will  bypass  this  option  and  not  try  to  open  and  read 
a  data  file.  After  receiving  the  input,  the  program  will  run.  It  will  echo 
the  input  data  to  the  screen,  with  several  PAUSES.  Press  the  ENTER  key  at 
each  PAUSE  to  continue  execution  of  the  program.  For  the  input  data  of 
Table  5,  it  turns  out  that  the  initial  guess  is  too  small  and  the  ranges  are 
small  also.  The  program  will  try  to  change  this  initial  guess  into  a  design 


28 


D:\>BUSHING 

ENTER  INPUT  FILE  NAME  PREFIX  FOR  A  .IN  ASCII  FILE. 

IT  MAY  INCLUDE  A  PATH  TO  ANOTHER  DIRECTORY, 

WITH  NO  MORE  THAN  26  CHARACTERS. 

IF  THE  FIRST  CONSTANT  IS  THE  NUMBER  OF  DATA  POINTS, 

AND  ITS  DESCRIPTION  INCLUDES  "DATA  POINTS." 

A  DATA  FILE  WITH  X,Y  DATA  PAIRS  AND  A  .DAT 
EXTENSION  SHOULD  ALSO  EXIST  IN  THE  SAME  DIRECTORY. 

THE  OUTPUT  FILE  WILL  HAVE  THE  SAME  PREFIX  WITH  A  .OUT  EXTENSION. 

NYLONl 


Seek  Screen  Image  at  Start  of  Program 
FIGURE  4 
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which  does  satisfy  the  design  constraints,  but  since  the  design  ranges  are 
small  also,  the  improveitient  steps  are  too  small  to  reach  the  acceptable  design 
space  in  the  twenty  steps  allowed  by  the  program. 

Table  6  is  a  listing  of  the  output  file  for  this  trial.  In  this  output, 
one  can  see  that  the  shaft  diameter  was  tripled  in  an  attempt  to  reach  the 
region  of  good  designs.  The  PV  limit  was  the  major  constraint,  with  values 
twenty  percent  higher  than  the  limit  for  the  revised  design  and  three  times 
the  limit  for  the  initial  design.  A  second  input  file,  NYL0N2.IN  can  now  be 
made  by  copying  the  first  and  modifying  the  design  parameter  initial  values. 

In  this  second  file  the  ranges  of  both  design  parameters  are  left  at  zero  to 
ten  mm  and  the  initial  values  are  increased  to  twenty  mm  for  the  length  and 
thirty  mm  for  the  diameter.  These  are  the  only  changes  from  the  first  input 
file.  Table  7  lists  the  new  input  file,  NYL0N2.IN. 

The  results  of  running  the  program  again  with  the  new  input  file  are 
listed  in  Table  8.  An  optimum  design  was  found  in  45  steps  with  a  length  of 
14.28  mm  and  a  diameter  of  20.41  mm.  The  design  has  a  friction  torque  of 
1.53  N-m  and  satisfies  all  three  constraints.  The  two  limiting  constraints 
which  are  just  satisfied  are  the  PV  factor  limit  and  the  length  to  diameter 
ratio.  These  results  are  consistent  with  the  graphical  results  of  Figures  2 
and  3. 

Note  that  the  input  ranges  of  the  design  variables  only  set  the 
sensitivity  of  the  search,  they  do  not  limit  the  design  to  have  diameters  less 
than  10  mm.  To  make  that  a  limit,  one  must  add  an  upper  limit  of  10  mm  to  the 
diameter  as  a  fourth  constraint.  This  would  prevent  the  optimizer  from 
finding  a  solution,  since  no  design  with  a  diameter  less  than  10  mm  can  have  a 
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Table  6 


First  Bushing  Output  Log  File 


RADIAL  NYLON  BUSHING 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 

FIXED  DESIGN  REQUIREMENTS: 

1  RADIAL  LOAD  750.00000  NEWTONS 

2  SHAFT  SPEED  40.00000  RPM 

3  FRICTION  COEFFICIENT  0.20000 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 

LOW  HIGH  INITIAL 


I  BUSHING  LENGTH 

0.00000 

10.00000 

5.00000  1 

2  BUSHING  DIAMETER 

0.00000 

10.00000 

5.00000  1 

THE  3  CONSTRAINT  LIMITS  ARE: 

TYPE 

I  CONTACT  PRESSURE 

14.00000 

MPa 

UPPER 

2  PV  FACTOR 

O.IIOOO 

MPa  -  m/s 

UPPER 

3  LENGTH  TO  DIAMETER 

0.70000 

RATIO 

UPPER 

MINIMIZE  THE  OBJECTIVE  FUNCTION. 

OBJ  =  FRICTION  TORQUE  IN  N  -  m  TIMES  I. 0000 
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Table  6  Continued 
First  Bushing  Output  Log  File 


OPTIMIZATION  COULD  NOT  BEGIN  -  BAD  INITIAL  VALUE. 


THE  INITIAL  VALUE  FOR  THE  PROBLEM  VIOLATED 
AT  LEAST  ONE  INEQUALITY  CONSTRAINT. 

THE  PROGRAM  COULD  NOT  FIND  ANOTHER  VECTOR  WHICH 
SATISFIED  ALL  THE  INEQUALITY  CONSTRAINTS. 


1  BUSHING  LENGTH 

2  BUSHING  DIAMETER 


THE  3  CONSTRAINT  VALUES  FOR 


1  CONTACT  PRESSURE  =  30.000 

2  PV  FACTOR  =  .31416 

3  LENGTH  TO  DIAMETER  =  I. 0000 

THE  3  CONSTRAINT  VALUES  FOR 


1  CONTACT  PRESSURE  =  3.6880 

2  PV  FACTOR  =  .13004 

3  LENGTH  TO  DIAMETER  =  .71744 


X  INITIAL 

X  MODIFIED 

5.00000 

12.07892  mm 

5.00000 

16.83625  mm 

INITIAL  ARE: 

LIMIT 

TYPE 

MPa 

14.000 

UPPER 

MPa  -  m/s 

.11000 

UPPER 

RATIO 

.70000 

UPPER 

MODIFIED  ARE: 

LIMIT 

TYPE 

MPa 

14.000 

UPPER 

MPa  -  m/s 

.11000 

UPPER 

RATIO 

.70000 

UPPER 
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Table  7 


Second  Bushing  Input  File 


RADIAL  NYLON  BUSHING 
3 

750.0 

RADIAL  LOAD 

NEWTONS 

40.0 

SHAFT  SPEED 
RPM 
0.2 

FRICTION  COEFFICIENT 
2 

0.0  10.0  20.0 
BUSHING  LENGTH 
mm 

0.0  10.0  30.0 

BUSHING  DIAMETER 
mm 
3 

UPPER  14.0 
CONTACT  PRESSURE 
MPa 

UPPER  0.11 
PV  FACTOR 
MPa  -  m/s 
UPPER  0.7 
LENGTH  TO  DIAMETER 
RATIO 
MIN 
1 

1.0 

FRICTION  TORQUE 
N  -  m 
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Table  8 


Second  Bushina  Output  Log  File 


RADIAL  NYLON  BUSHING 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 

FIXED  DESIGN  REQUIREMENTS: 

1  RADIAL  LOAD  750.00000  NEWTONS 

2  SHAFT  SPEED  40.00000  RPM 

3  FRICTION  COEFFICIENT  0.20000 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 


LOW 

HIGH 

INITIAL 

1  BUSHING  LENGTH 

0.00000 

10.00000 

20.00000  mm 

2  BUSHING  DIAMETER 

0.00000 

10.00000 

30.00000  mm 

THE  3  CONSTRAINT  LIMITS  ARE: 

TYPE 

1  CONTACT  PRESSURE 

14.00000 

MPa 

UPPER 

2  PV  FACTOR 

0.11000 

MPa  -  m/s 

UPPER 

3  LENGTH  TO  DIAMETER 

0.70000 

RATIO 

UPPER 

MINIMIZE  THE  OBJECTIVE  FUNCTION. 

OBJ  =  FRICTION  TORQUE 

IN  N 

-  m 

TIMES  1.1 

OPTIMIZATION  SUCCESSFUL  IN  45  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 

X(I) 

1  BUSHING  LENGTH  14.28342  mm 

2  BUSHING  DIAMETER  20.40649  mm 

THE  MINIMUM  OBJECTIVE  FUNCTION  =  1.53049  ,  ITS  COMPONENTS  ARE: 

1  FRICTION  TORQUE  =  1.5305  N  -  m  TIMES 

THE  LAST  CHANGE  IN  THE  OBJECTIVE  FUNCTION  =  -0.677109E-04 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.390625E-03 


.0000 
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Table  8  Continued 


THE  3  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1  CONTACT  PRESSURE 

=  2.5735 

MPa 

14.000 

UPPER 

2  PV  FACTOR 

=  .10999 

MPa  -  m/s 

.11000 

UPPER 

3  LENGTH  TO  DIAMETER 

=  .69989 

RATIO 

.70000 

UPPER 

DESIGN  CHECK 

X(I) 

1  BUSHING  LENGTH  15.00000  mm 

2  BUSHING  DIAMETER  21.00000  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  1.57500  ,  ITS  COMPONENTS  ARE: 

1  FRICTION  TORQUE  =  1.5750  N  -  m  TIMES  1.0000 


THE  3  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1  CONTACT  PRESSURE 

=  2.3810 

MPa 

14.000 

UPPER 

2  PV  FACTOR 

=  .10472 

MPa  -  m/s 

.11000 

UPPER 

3  LENGTH  TO  DIAMETER 

=  .71429 

RATIO 

.70000 

UPPER 

DESIGN  CHECK 

X(I) 

1  BUSHING  LENGTH  15.00000  mm 

2  BUSHING  DIAMETER  22.00000  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  1.65000  ,  ITS  COMPONENTS  ARE: 

1  FRICTION  TORQUE  =  1.6500  N  -  m  TIMES  1.0000 


THE  3  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1  CONTACT  PRESSURE 

=  2.2727 

MPa 

14.000 

UPPER 

2  PV  FACTOR 

=  .10472 

MPa  -  m/s 

.11000 

UPPER 

3  LENGTH  TO  DIAMETER 

=  .68182 

RATIO 

.70000 

UPPER 
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PV  factor  less  than  0.11  MPa  -  m/s  and  a  P  less  than  0.7.  Repeated  trials 
such  as  the  first  one  would  tell  us  that. 

However,  another  constraint  on  the  solution  was  that  the  diameter  and 
length  be  in  whole  mm’s.  This  can  be  obtained  with  the  design  check  provision 
of  the  program,  which  is  shown  in  Table  8.  Once  the  numerical  optimum  has 
been  found,  the  program  lists  the  number  of  optimizing  steps  followed  by:  the 
found  design  parameter  values  for  bushing  length  and  bushing  diameter,  the 
objective  function  value  and  the  three  constrained  variable  values  and  limits. 
Then  the  program  re-lists  the  design  variables  with  their  found  values  and 
offers  the  user  the  option  to  change  them  for  a  design  check.  Figure  5  shows 
this  interaction.  The  user  responded  with  a  ’Y’  to  the  question  on  trying 
another  design  and  entered  the  two  values  of  ’15’  and  ’21’  for  the  bushing 
length  and  diameter.  The  program  then  printed  the  results  to  the  screen  and 
added  them  to  the  output  file  as  shown  in  Table  8.  This  option  is  offered  to 
the  user  at  the  end  of  each  analysis  until  the  response  to  the  first  question 
is  ’N’  which  tells  the  program  to  close  the  output  file  and  stop  the  program. 

Increasing  the  length  to  15  mm  and  the  diameter  to  21  mm  increases  the 
friction  torque  slightly  to  1.575  N-m  but  still  does  not  satisfy  the  length  to 
diameter  constraint.  A  second  trial  with  a  15  mm  length  and  a  22  mm  diameter 
has  a  friction  torque  of  1.65  mm  and  satisfies  all  three  modeled  constraints 
and  the  additional  requirement  of  standard  sizes.  This  trial  is  the  optimal 
design  and  is  shown  in  Figure  6. 
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FOR  THE  LAST  DESIGN,  THE  DESIGN  VARIABLE  VECTOR  IS: 

X{I) 

1  BUSHING  LENGTH  14.28342  mm 

2  BUSHING  DIAMETER  20.40649  mm 


DO  YOU  WISH  TO  MODIFY  THE  DESIGN  AND  CHECK  IT  ?  (Y/N) 


ENTER  NEW  VALUES  FOR  THE  DESIGN  VARIABLE  VECTOR,  X. 
ENTER  A  COMMA  TO  KEEP  THIS  VALUE. 


OLD  BUSHING  LENGTH  =  14.28342  mm 

NEW  BUSHING  LENGTH  =  ? 

15 

OLD  BUSHING  DIAMETER  =  20.40649  mm 

NEW  BUSHING  DIAMETER  =  ? 


Seek  Design  Verification  Request  Screen 
FIGURE  5 
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1  5  mm 


FIGURE  6 


SPRING  OPTIMIZATION 


In  this  example,  consider  the  design  of  a  steel  helical  coil  compression 
spring  to  support  a  varying  load  with  a  specified  spring  rate.  The  applied 
load  varies  between  minimum  and  maximum  values  in  service  and  may  be  larger  if 
the  spring  is  compressed  solid  at  assembly.  For  a  given  spring  material,  four 
geometric  parameters  define  the  spring:  1)  the  wire  diameter,  d^;  2)  the  mean 
coil  diameter,  D;  3)  the  number  of  active  coils,  N^^;  and  4)  the  height  of  the 
spring  when  unloaded,  h^. 

In  this  design  problem,  one  more  requirement  can  be  placed  on  the 
performance  of  the  spring:  it  could  have  a  specified  outside  diameter,  or  work 
over  a  rod  of  a  given  diameter  or  it  could  have  a  required  height  under  load. 
Instead,  we  will  let  the  optimizer  find  a  spring  with  a  minimized  property 
which  can  support  the  specified  loads  with  the  given  spring  rate.  Three 
separate  objective  functions  will  be  minimized:  1)  spring  weight,  2)  spring 
height,  and  3)  spring  coil  volume.  The  resulting  designs  will  satisfy  the 
load  and  deflection  requirements,  but  they  will  be  different. 

For  the  spring  design  to  have  some  practicality,  the  wire  size  should  be 
selected  from  a  finite  list  of  available  diameters  and  the  mean  coil  diameter 
should  also  be  a  standard  size. 

Theory 

A  spring’s  performance  is  modeled  by  its  strength  and  deflection.  A 
helical  coil  spring  supports  its  axial  load  as  a  torsional  shear  stress  in  the 
wire  with  a  small  additional  direct  shear  stress.  Figure  7  shows  the  axial 
load  and  the  internal  wire  torque  and  shear  which  support  it.  Due  to  the 
curvature  of  the  wire,  Wahl  determined  a  stress  concentration  factor, 
which  compares  the  maximum  shear  stress  in  the  wire  to  the  nominal  torsional 
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SHEAR  STRESSES 

HELICAL  COIL  SHEAR  LOADS  AND  STRESSES 

FIGURE? 
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stress  caused  by  the  couple  of  the  axial  load  through  the  center  of  the  spring 
and  the  supporting  shear  force  through  the  center  of  the  wire  [10,15].  The 
shear  stress  and  stress  concentration  factor  are: 


T  = 


Ku,  T,  D 
w  q 

2.0  J 


8.0  K,,  F  D 
w 


n  d 


w 


and 


K. 


4.0  C  -  1.0  0.615 

-  +  - 


4.0  C  -  4.0 

where  the  spring  index,  C  is: 

D 

C  =  - 


w 


8.0  K,  F  C 
w 


(8) 


(9) 


(10) 


The  deflection  of  a  helical  coil  spring  at  any  load,  F,  is  a  result  of 
the  twist  in  the  wire  due  to  the  applied  torque,  T^: 


D  D  L 

6  =  —e  = - 9— 


2  J  G 


(11) 


or 


or 


6  = 


D/2  (  F  D/2  )  (  77  D  ) 
(  n  d^VsZ  )  G 


8  F  r  N. 


8  F  D-"  N. 
.  4  g  ‘ 


6  = 


(12) 


(13) 


The  stiffness  of  the  spring  or  its  rate  is  thus: 


8.0  r  N. 


(14) 
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These  basic  equations  for  stress  (8),  deflection  (13)  and  stiffness  (14) 
can  be  applied  to  the  specified  loads  to  model  the  performance  of  the  spring. 
To  avoid  yielding,  the  maximum  stress  in  the  spring  should  be  less  than  the 
shear  yield  strength  of  the  wire  divided  by  the  design  stress,  or  the  static 
design  factor,  N^,  should  be  greater  than  the  desired  design  factor,  where: 

N  .  (15) 

^  '’sol 

and  the  maximum  stress  in  the  spring,  is  found  by  using  the  solid  height 

force,  ,  in  equation  (8).  The  shear  yield  strength  of  the  wire  can  be 
estimated  as  eighty  percent  of  the  ultimate  shear  strength  of  the  wire.  As 
spring  wire  is  drawn  to  a  smaller  diameter,  its  ultimate  strength  increases. 
Using  the  octahedral  shear  factor  of  0.577,  the  ultimate  shear  strength  can  be 
estimated  as  a  function  of  the  drawn  wire  size  as  [15]: 

S  =  0.577  d^^  (16) 

su  uc  w 

In  addition  to  having  adequate  strength  to  avoid  yielding  when 
compressed  solid,  springs  should  be  able  to  support  an  infinite  number  of 
loading  cycles  without  experiencing  a  fatigue  failure.  Figure  8  is  a 
Soderberg  diagram  of  alternating  shear  stress  versus  mean  shear  stress,  which 
shows  the  reduction  in  alternating  strength  with  increasing  mean  stress,  for  a 
helical  coil  compression  spring.  Since  the  stress  concentration  affects  the 
alternating  stress  directly  and  does  not  affect  the  mean  stress  influence,  the 
Wahl  factor  is  applied  only  to  the  alternating  shear  stress.  An  algebraic 
trick,  which  allows  all  stresses  to  be  calculated  with  equation  (8),  is  to 
multiply  the  stress  concentration  factor  by  the  ultimate  strength  to  divide  it 
out  of  the  mean  stress. 
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ALTERNATING  SHEAR  STRESS  (MPa) 


In  Figure  8,  the  negative  sloped  lines  show  the  Goodman  criteria  for 

fatigue  strength  with  and  without  the  design  factor,  N^,  and  the  positive 

sloped  line  shows  the  load  ratio  of  alternating  stress,  r^,  to  mean  stress, 

T  for  this  application.  The  design  factor  equation  for  this  criteria  is 
m’ 

[15]: 


N, 


*^w  ^su 


(17) 


se 


or 


1.0 


(18) 


Kw'su 


se 


The  mean  stress  is  calculated  with  equation  (8)  using  the  average  applied 
load; 


F  +  F  . 
max  mm 


(19) 


2.0 

And  the  alternating  stress  is  calculated  using  the  alternating  applied  load: 

(20) 


max 


^min 


2.0 


Unlike  the  ultimate  wire  strength,  the  fatigue  strength  of  the  steel  wire, 

S  ,  is  constant  at  these  high  strengths, 
se 

Figure  9  is  a  plot  of  the  force  on  the  spring  versus  the  height  of  the 
spring.  Near  the  solid  height,  the  axial  load  on  the  spring  increases  rapidly 
as  the  coils  close.  In  service,  the  spring  operates  between  the  working 
heights  of  h^^.^  and  with  loads  of  F^^^  and  The  unloaded  height  of 

the  spring  is  the  sum  of  the  deflection  to  solid  height  and  the  solid  height. 
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AXIAL  LOAD -F  (N) 


^solid  ^min  ^max 

SPRING  HEIGHT  -  h  (mm) 

SPRING  FORCE  VERSUS  POSITION 
FIGURE  9 
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sol 


+  1.01  (  Ng  +  Ng  ) 


(21) 


Coil  weight  is  a  direct  function  of  the  volume  of  wire  in  the  spring: 


n 


w 


d„M  N,  +  )  >7  D 


(22) 


with: 


Wt  =  wV„ 


(23) 


The  volume  of  the  coil  is  the  area  of  the  outside  diameter’s  circular  disc 
times  the  spring  height: 


n  2 

V  =  —  OD^  hr 

coil  4  f 


(24) 


Programming 

Table  9  summarizes  the  design  problem  in  terms  of  the  constants  which 
define  the  problem,  the  design  parameters  which  are  to  be  found,  the  equality 
and  inequality  constraints  on  the  design  and  the  three  separate  objective 
functions  which  will  be  sought.  In  the  initial  formulation,  there  are  two 
equality  constraints  on  the  stiffness  and  the  force  at  solid  height,  which  can 
be  used  to  reduce  the  number  of  independent  design  parameters  from  four  to 
two.  There  are  four  active  inequality  constraints:  1)  that  the  number  of 
active  coils  be  positive,  2)  that  the  fatigue  design  factor  be  greater  than 
the  desired  design  factor,  3)  that  the  static  design  factor  at  assembly  be 
greater  than  the  desired  design  factor,  and  4)  that  the  spring  index  be 
greater  than  two,  so  a  coil  can  be  manufactured. 
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Table  9 


Initial  Soring  Optimization  Parameters 


Constants 


Design  Parameters 


Constraints 
Equality  Inequality 


Objective  Function 


min 

^max 

^solid 


k 

N, 


"des 

^se 

^su 


G 

w 


w 


D 

N. 


solid 


N,  >  0.0 

a 


Nf  >  "des 


Ns  >  "des 


C  >  2.0 


^Mmin 


or 


'  t^min 


or 


'  'min 


Table  10 

Revised  Soring  Optimization  Parameters 


Constants 

Design  Parameters 

Inequality  Constraints 

Objective  Function 

^min 

N,  >  0.0 
a 

^•’f^min 

*^max 

D 

>  Ndes 

or 

*^sol  id 

“s  >  "des 

'  t'min 

k 

C  >  2.0 

or 

Ne 

>0.0 

max 

(Vol)„d„ 

'  'min 

s 

OD  >  0.0 

se 

s 

Vw  >  0-0 

uc 

w 

a 

G 

w 
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From  the  stiffness  equation  (14),  one  can  relate  the  number  of  active 
coils  to  the  wire  and  mean  coil  diameters,  the  stiffness  and  the  shear  modulus 
of  the  material: 


(25) 


Equation  (21)  relates  the  spring  height  to:  the  maximum  force,  the 

spring  rate;  the  wire  diameter;  and  the  numbers  of  active  and  end  coils. 

These  two  equations  convert  two  design  parameters  from  independent  parameters 
to  dependent  parameters  and  simplify  the  optimization.  Table  10  is  a  second 
pass  at  formulating  the  optimization  problem  in  this  simpler  form.  Three 
inactive  constraints  have  been  added  to  the  inequality  constraint  list  to  make 
the  optimizer  report  these  properties.  All  properties  are  constrained  to  be 
greater  than  zero,  which  they  will  be  for  all  designs.  The  watched  properties 
are:  1)  maximum  static  stress,  2)  outside  diameter,  OD;  and  3)  wire 

volume,  V  .  The  final  differences  between  the  two  formulations  are  in  the 
constant  list:  1)  the  elimination  of  the  design  factor  from  the  constant  list 
since  it  is  used  in  the  constraint  limit  list  and,  2)  the  replacement  of  the 
ultimate  shear  strength  of  the  wire  by  the  ultimate  tensile  strength  constant 
and  the  wire  power  to  enable  the  program  to  vary  the  strengths  with  wire  size. 

Figure  10  is  a  plot  of  the  design  space  for  this  reformulated 
optimization.  The  graph  plots  the  two  independent  design  parameters:  the  wire 
diameter,  d  ,  and  the  mean  coil  diameter,  D,  versus  each  other.  On  the  graph 
are  drawn  constraint  lines  for  the  four  active  constraints.  The  number  of 
active  coils  constraint,  N  ,  is  drawn  for  a  limit  of  0.1  coils  to  show  its 

a 

shape.  A  limit  of  zero  coils  lies  on  the  x  axis.  Designs  which  satisfy  all 
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HELICAL  COIL  SPRING  DESIGN  SPACE 

FIGURE  10 


constraints  can  be  found  in  the  region  labeled  acceptable  designs  which  is 
between  the  spring  index  limit,  labeled  C,  and  the  solid  stress  limit,  labeled 

SOLID. 

Each  of  the  three  objective  functions  are  complex  functions  of  these 
parameters,  so  a  single  plot  of  the  objective  function  versus  mean  coil 
diameter  is  not  possible.  Contour  plots  for  each  objective  function  would 
have  to  be  drawn  on  the  design  space  to  obtain  a  graphical  solution  to  the 

problem. 

Two  analysis  subroutines  must  now  be  written  to  operate  on  the  constants 
and  design  parameters  of  Table  10  and  determine  the  constraint  values  and  the 
objective  function  values  listed.  Subroutine  BOUNDS,  which  is  listed  in 
Table  11,  takes  the  constant  array  and  the  design  parameter  array  from  the 
calling  list  and  determines  the  constraint  values.  As  with  the  bushing 
example,  this  is  done  in  three  steps  to  clarify  the  calculations: 

1)  conversion  of  input  arrays  to  variables,  2)  calculation  of  the  properties 
and  3)  transferring  the  property  values  to  the  constraint  property  array.  The 
analysis  of  equations  (8)  through  (25)  is  used  in  the  subroutine.  Subroutine 
VALUES,  which  is  listed  in  Table  12,  performs  a  similar  determination  of  the 
three  objective  function  values.  Once  written,  these  two  subroutines  must  be 
compiled  and  linked  to  program  SEEK  to  generate  an  executable  program  to 
perform  the  optimization. 

Numerical  Example 

For  the  example,  consider  the  design  of  a  spring  to  have  a  spring  rate 
of  4  kN/m.  In  use,  the  spring  will  see  a  load  which  varies  from  a  minimum  of 
24  N  to  a  maximum  of  120  N.  The  spring  should  not  go  solid  until  the  applied 
load  reaches  a  value  of  150  N.  The  spring  wire  is  to  be  selected  from  stock 
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Table  11 


Spring  Constraint  Evaluation  Subroutine  Bounds 


SUBROUTINE  BOUNDS ( CONST , NCO , X , NX , VCSTR , NCS ) 

BOUNDS  DETERMINES  THE  PRESENT  CONSTRAINT 
FUNCTION  VALUES 

FOR  A  HELICAL  COIL  SPRING  DESIGN 
PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NCS  -  NUMBER  OF  INEQUALITY  CONSTRAINTS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  VARIABLES 

VCSTR  -  PRESENT  CONSTRAINT  VALUES 

X  -  PRESENT  DESIGN  VARIABLE  VALUES 

ALL  VALUES  ARE  IN  PROBLEM  UNITS 


CONST(l)  =  Fmin 
CONST(2)  =  Fmax 
CONST(3)  =  Fsol 
CONST(4)  =  k 
CONST(5)  =  Ne 
C0NST(6)  =  Sse 
CONST(7)  =  Sue 
C0NST(8)  =  a 
C0NST(9)  =  G 
CONST(IO)  =  w 


MINIMUM  FORCE  (NEWTONS) 

MAXIMUM  FORCE  (NEWTONS) 

MINIMUM  SOLID  FORCE  (NEWTONS) 
SPRING  RATE  (kN/m) 

NUMBER  OF  DEAD  COILS 
SHEAR  FATIGUE  STRENGTH  (MPa) 
TENSILE  STRENGTH  COEFICIENT  (MPa) 
TENSILE  STRENGTH  WIRE  POWER 
SHEAR  MODULUS  (MPa) 

WEIGHT  DENSITY  (kN/m**3) 


X(l)  =  dw  -  WIRE  DIAMETER  (mm) 

X(2)  =  D  -  MEAN  COIL  DIAMETER  (mm) 


VCSTR(l)  =  Na 
VCSTR(2)  =  Nf 
VCSTR(3)  =  Ns 
VCSTR(4)  =  C 
VCSTR(5)  =  Tmax 
VCSTR(6)  =  OD 
VCSTR(7)  =  Vw 


NUMBER  OF  ACTIVE  COILS 
FATIGUE  DESIGN  FACTOR 
STATIC  DESIGN  FACTOR 
SPRING  INDEX 

MAXIMUM  SHEAR  STRESS  (MPa) 
OUTSIDE  DIAMETER  (mm) 
SPRING  WIRE  VOLUME  (mm**3) 
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Table  11  Continued 

Soring  Constraint  Evaluation  Subroutine  Bounds 


DIMENSION  CONST(NCO) ,X(NX) ,VCSTR(NCS) 

PI  =  3.14159265 
FMIN  =  CONST(l) 

FMAX  =  C0NST(2) 

FSOL  =  C0NST(3) 

RATE  =  CONST(4) 

ZNE  =  CONST(5) 

SSE  =  CONST(6) 

CSUT  =  C0NST(7) 

ASUT  =  C0NST(8) 

G  =  C0NST(9) 

DW  =  X(l) 

D  =  X(2) 

C  =  D/DW 

ZNA  =  DW*G/(8.0*C*C*C*RATE) 

SSU  =  0.577  *  CSUT  *  DW**ASUT 
SSY  =  0.8*SSU 
FA  =  (FMAX  -  FMlN)/2.0 
FM  =  (FMAX  +  FMIN)/2.0 

SKW  =  (4.0*C  -  1.0)/(4.0*C  -  4.0)  +  0.615/C 
TA  =  (8.0*SKW*FA*C/(PI*DW*DW)) 

TM  =  (8.0*FM*C/(PI*DW*DW)) 

ZNF  =  1.0/(TM/SSU  +  TA/SSE) 

TMAX  =  (8.0*SKW*FSOL*C/(PI*DW*DW)) 

ZNS  =  SSY/TMAX 

OD  =  D  +  DW 

ZNT  =  ZNA  +  ZNE 

VW  =  0.25*PI*DW*DW*ZNT*PI*D 

VCSTR(l)  =  ZNA 

VCSTR(2)  =  ZNF 

VCSTR(3)  =  ZNS 

VCSTR(4)  =  C 

VCSTR(5)  =  TMAX 

VCSTR(6)  =  OD 

VCSTR(7)  =  VW 

RETURN 

END 
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Table  12 


Soring  Objective  Function  Evaluation  Subroutine  Values 


SUBROUTINE  VALUES(CONST,NCO,X, NX, OBJECT, NOB) 

VALUES  DETERMINES  THE  PRESENT  DESIGN 
OBJECTIVE  FUNCTION  VALUES 

FOR  A  CANTILEVER  BEAM  DESIGN  EXAMPLE 


PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NOB  -  NUMBER  OF  OBJECTIVE  FUNCTION  TERMS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  VARIABLES 

OBJECT  -  PRESENT  OBJECTIVE  FUNCTION  VALUES 

X  -  PRESENT  DESIGN  VARIABLE  VALUES 


ALL  VALUES  ARE  IN  PROBLEM  UNITS 


CONST(l)  =  Fmin 
CONST(2)  =  Fmax 
CONST(3)  =  Fsol 
CONST(4)  =  k 
CONST(5)  =  Ne 
C0NST(6)  =  Sse 
CONST(7)  =  Sue 
C0NST(8)  =  a 
C0NST(9)  =  G 
CONST(IO)  =  w 


MINIMUM  FORCE  (NEWTONS) 

MAXIMUM  FORCE  (NEWTONS) 

MINIMUM  SOLID  FORCE  (NEWTONS) 
SPRING  RATE  (kN/m) 

NUMBER  OF  DEAD  COILS 
SHEAR  FATIGUE  STRENGTH  (MPa) 
TENSILE  STRENGTH  COEFICIENT  (MPa) 
TENSILE  STRENGTH  WIRE  POWER 
SHEAR  MODULUS  (MPa) 

WEIGHT  DENSITY  (kN/m**3) 


X(l)  =  dw 
X(2)  =  D 


WIRE  DIAMETER  (mm) 

MEAN  COIL  DIAMETER  (mm) 


OBJECT(l)  =  Wt  -  SPRING  WEIGHT  (NEWTONS) 

OBJECT(2)  =  hf  -  SPRING  HEIGHT  (mm) 

OBJECT(3)  =  Vc  -  SPRING  CYLINDER  VOLUME  (mm**3) 


DIMENSION  CONST(NCO) ,X(NX) ,OBJECT(NOB) 
PI  =  3.14159265 
FSOL  =  CONST(3) 

RATE  =  CONST (4) 

ZNE  =  CONST(5) 

G  =  CONST(9) 

W  =  C0NST(10)/1000000.0 
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Table  12  Continued 

Soring  Objective  Function  Evaluation  Subroutine  Values 


DSOL  =  FSOL/RATE 
DW  =  X(l) 

D  =  X(2) 

C  =  D/DW 
OD  =  D  +  DW 

ZNA  =  DW*G/(8.0*C*C*C*RATE) 

ZNT  =  ZNA  +  ZNE 

VW  =  0.25*PI*DW*DW*ZNT*PI*D 

Wj  =  vw*w 

HE  =  DSOL  +  DW*ZNT*1.01 
VC  =  0.25*PI*0D*0D*HF 
OBJECT(l)  =  WT 
OBJECT(2)  =  HF 
OBJECT(3)  =  VC 
RETURN 
END 
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sizes  which  have  whole  and  half  mm  values.  The  mean  coil  diameter  should  also 
have  values  with  whole  or  half  mm  precision.  The  springs  are  to  be  made  of 
hard  drawn  spring  wire  with  an  ultimate  strength  constant  of  1510  MPa,  a  wire 
diameter  strength  variation  exponent  of  -0.201,  and  a  shear  fatigue  strength 
of  465  MPa.  The  acceptable  design  factor  is  1.5  and  the  spring  ends  are  to  be 
squared  and  ground  with  one  inactive  coil  at  each  end.  The  shear  modulus  of 

3 

steel  is  79,000  MPa  and  its  weight  density  is  76.5  kN/m  . 

Table  13  is  a  listing  of  the  input  file  for  the  minimum  weight  design 
option.  The  file  begins  with  the  title  for  the  output  file.  The  second  line 
has  the  number  of  constants  to  follow  -  ten.  The  next  thirty  lines  contain 
these  ten  constants  in  the  order  listed  in  Table  10  with  their  descriptions 
and  units.  Following  the  constants  is  a  line  with  the  number  two  which 
indicates  that  two  independent  design  parameters  will  follow.  These  two 
parameters  are  the  wire  diameter  and  the  mean  coil  diameter.  Low,  high  and 
initial  estimates  are  selected  as  1.0,  15.0  and  5.0  for  the  wire  size  and 
25.0,  500.0  and  100.0  for  the  mean  coil  diameter.  The  next  line  has  the 
single  value  of  seven  for  the  seven  design  constraints  listed  in  Table  10. 

All  seven  constraints  happen  to  be  lower.  Following  the  constraints  is  a  line 
with  the  letters  ’MIN’  to  select  minimization  as  the  direction  of 
optimization,  a  line  with  the  value  of  three  to  indicate  that  there  are  three 
terms  in  the  objective  function.  The  last  nine  lines  are  the  weighting 
coefficients,  names  and  units  for  the  three  objective  function  terms  of 
weight,  height  and  volume.  The  weighting  coefficient  of  the  first  is  one  and 
the  last  two  are  zero.  By  changing  which  term  has  the  unit  coefficient,  one 
can  change  the  optimization  goal  without  changing  the  program. 
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Table  13 


Soring  Design  Input  File  For  Minimum  Weight 


HELICAL  COIL  SPRING  -  MINIMUM  WEIGHT 
10 

24.0 

MINIMUM  FORCE 
NEWTONS 
120.0 

MAXIMUM  FORCE 

NEWTONS 

150.0 

MINIMUM  FORCE  WHEN  SOLID 

NEWTONS 

4.0 

SPRING  RATE 

kN/m 

2 

END  COIL  NO. 

465.0 

SHEAR  FATIGUE  STRENGTH 

MPa 

1510.0 

TENSILE  STRENGTH  CONSTANT 
MPa 
-0.201 

TENSILE  STRENGTH  WIRE  POWER 

79000.0 
SHEAR  MODULUS 
MPa 
76.5 

WEIGHT  DENSITY 
kN/ni**3 
2 

1.0  15.0  5.0 

WIRE  DIAMETER 
mm 

25.0  500.0  100.0 

MEAN  COIL  DIAMETER 
mm 
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Table  13  Continued 


Spring  Design  Input  File  For  Minimum  Weight 


7 

LOWER  0.0 
ACTIVE  COIL 
NO. 

LOWER  1.5 

FATIGUE  DESIGN  FACTOR 
LOWER  1.5 

STATIC  DESIGN  FACTOR 

LOWER  2.0 
SPRING  INDEX 

LOWER  0.0 

MAXIMUM  SHEAR  STRESS 
MPa 

LOWER  0.0 

OUTSIDE  COIL  DIAMETER 
mm 

LOWER  0.0 

SPRING  WIRE  VOLUME 

mm**3 

MIN 

3 

1.0 

WEIGHT 

NEWTONS 

0.0 

HEIGHT 

mm 

0.0 

SPRING  COIL  VOLUME 
mm**3 
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This  ASCII  file  can  then  be  saved  with  a  name  such  as  ’WEIGHT. IN’  and 

used  with  the  spring  optimization  program.  Table  14  is  the  output  file  which 

resulted  from  running  the  spring  optimization  program  with  this  file.  The  ten 

constants,  two  design  parameters,  seven  constraints  and  three  objective 

function  term  are  all  listed  with  their  values,  names,  units,  limit 

directions,  and  weighting  coefficients  at  the  start  of  the  file.  The 

optimization  reached  a  solution  in  24  steps  with  a  wire  diameter  of  slightly 

more  than  3  mm  and  a  mean  coil  diameter  of  23.13  mm.  The  spring  weight  was 

0.79  Newtons  and  the  spring  had  17.5  active  coils  with  a  height  of  97  mm  a 

3 

spring  index  of  7.5  and  a  spring  coil  volume  of  52,000  mm  .  The  static 
overload  stress  was  the  limiting  factor  in  the  design  with  a  static  design 
factor  of  1.5  at  the  limit  and  a  maximum  shear  stress  of  370  MPa. 

Following  this  output  is  a  design  check  with  a  3.5  mm  wire  diameter  and 
a  mean  coil  diameter  of  35  mm.  With  the  larger  standard  wire  size,  the  larger 
coil  diameter  keeps  the  spring  weight  down  by  reducing  the  number  of  active 
coils  needed  to  obtain  the  spring  rate  without  lowering  the  design  factors. 

The  weight  increased  to  0.86  Newtons,  the  number  of  active  coils  dropped  to 
8.64  and  the  spring  index  increased  to  10.  In  addition,  the  height  dropped  to 
75  mm  and  the  spring  coil  volume  increased  to  87,500  mm^.  Although  the  spring 
is  slightly  heavier  than  the  initial  optimum,  it  has  standard  wire  and  coil 
dimensions,  so  it  is  a  practical  optimal  solution  to  the  posed  problem. 

Table  15  is  the  output  file  produced  by  running  the  spring  optimization 
program  with  a  second  input  file  which  has  two  weighting  coefficients  changed 
and  the  problem  title  changed  in  line  one.  The  two  numerical  changes  were  to 
the  first  and  second  weighting  coefficient  values  to  switch  the  object  of 
minimization  from  weight  to  height.  As  shown  in  Table  15,  this  optimum  is 
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Table  14 


Soring  Design  Output  File  For  Minimum  Weight 


HELICAL  COIL  SPRING  -  MINIMUM  WEIGHT 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 


FIXED  DESIGN  REQUIREMENTS: 


1 

MINIMUM  FORCE 

24.00000 

NEWTONS 

2 

MAXIMUM  FORCE 

120.00000 

NEWTONS 

3 

MINIMUM  FORCE  WHEN  SOLID 

150.00000 

NEWTONS 

4 

SPRING  RATE 

4.00000 

kN/m 

5 

END  COIL  NO. 

2.00000 

6 

SHEAR  FATIGUE  STRENGTH 

465.00000 

MPa 

7 

TENSILE  STRENGTH  CONSTANT 

1510.00000 

MPa 

8 

TENSILE  STRENGTH  WIRE  POWER 

-0.20100 

9 

SHEAR  MODULUS 

79000.00000 

MPa 

10 

WEIGHT  DENSITY 

76.50000 

kN/m**3 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 

LOW  HIGH  INITIAL 

1  WIRE  DIAMETER  1.00000  15.00000  5.00000  mm 

2  MEAN  COIL  DIAMETER  25.00000  500.00000  100.00000  mm 

THE  7  CONSTRAINT  LIMITS  ARE: 

TYPE 


1 

ACTIVE  COIL 

0.00000 

NO. 

LOWER 

2 

FATIGUE  DESIGN  FACTOR 

1.50000 

LOWER 

3 

STATIC  DESIGN  FACTOR 

1.50000 

LOWER 

4 

SPRING  INDEX 

2.00000 

LOWER 

5 

MAXIMUM  SHEAR  STRESS 

0.00000 

MPa 

LOWER 

6 

OUTSIDE  COIL  DIAMETER 

0.00000 

mm 

LOWER 

7 

SPRING  WIRE  VOLUME 

0.00000 

mm**3 

LOWER 

MINIMIZE  THE  OBJECTIVE  FUNCTION,  WHICH 

HAS  3  TERMS. 

OBJ  =  THE  LINEAR  SUM  OF: 

1 

WEIGHT  IN 

NEWTONS 

TIMES 

1.0000 

2 

HEIGHT  IN 

mm 

TIMES 

0.0000 

3 

SPRING  COIL  VOLUME  IN 

mm**3 

TIMES 

0.0000 
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Table  14  Continued 

Soring  Design  Output  File  For  Minimum  Weight 


OPTIMIZATION  SUCCESSFUL  IN  24  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 


X(I) 


1  WIRE  DIAMETER  3.05491  mm 

2  MEAN  COIL  DIAMETER  23.12969  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  0.789480 


ITS  COMPONENTS  ARE: 


1  WEIGHT 

2  HEIGHT 

3  SPRING  COIL  VOLUME 


=0.78948 
=  97.285 
=  52388. 


NEWTONS 

mm 

mm**3 


TIMES 

TIMES 

TIMES 


1.0000 

0.0000 

0.0000 


THE  LAST  CHANGE  IN  THE  OBJECTIVE  FUNCTION  =  0.596046E-07 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.390625E-03 


THE  7  CONSTRAINT  VALUES  ARE: 


LIMIT  TYPE 


1  ACTIVE  COIL 

=  17.525 

NO. 

.00000 

2  FATIGUE  DESIGN  FACTOR 

=  2.1345 

1.5000 

3  STATIC  DESIGN  FACTOR 

=  1.5035 

1.5000 

4  SPRING  INDEX 

=  7.5462 

2.0000 

5  MAXIMUM  SHEAR  STRESS 

=  370.49 

MPa 

.00000 

6  OUTSIDE  COIL  DIAMETER 

=  26.070 

mm 

.00000 

7  SPRING  WIRE  VOLUME 

=  10320. 

mm**3 

.00000 

LOWER 

LOWER 

LOWER 

LOWER 

LOWER 

LOWER 

LOWER 
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Table  14  Continued 

Soring  Design  Output  File  For  Minimum  Weight 


DESIGN  CHECK 

X(I) 

1  WIRE  DIAMETER  3.50000  mm 

2  MEAN  COIL  DIAMETER  35.00000  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  0.861137  ,  ITS  COMPONENTS  ARE: 


1  WEIGHT 

=0.86114 

NEWTONS 

TIMES 

1.0000 

2  HEIGHT 

=  75.115 

mm 

TIMES 

0.0000 

3  SPRING  COIL  VOLUME 

=  87445. 

mm**3 

TIMES 

0.0000 

THE  7  CONSTRAINT 

VALUES  ARE: 

LIMIT 

TYPE 

1  ACTIVE  COIL 

=  8.6406 

NO. 

.00000 

LOWER 

2  FATIGUE  DESIGN  FACTOR 

=  2.1430 

1.5000 

LOWER 

3  STATIC  DESIGN  FACTOR 

=  1.5179 

1.5000 

LOWER 

4  SPRING  INDEX 

=  10.000 

2.0000 

LOWER 

5  MAXIMUM  SHEAR  STRESS 

=  356.97 

MPa 

.00000 

LOWER 

6  OUTSIDE  COIL  DIAMETER 

=  38.500 

mm 

.00000 

LOWER 

7  SPRING  WIRE  VOLUME 

=  11257. 

mm**3 

.00000 

LOWER 
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Table  15 


Spring  Design  Output  File  For  Minimum  Height 


HELICAL  COIL  SPRING  -  MINIMUM  HEIGHT 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 


FIXED  DESIGN  REQUIREMENTS: 


1 

MINIMUM  FORCE 

24.00000 

NEWTONS 

2 

MAXIMUM  FORCE 

120.00000 

NEWTONS 

3 

MINIMUM  FORCE  WHEN  SOLID 

150.00000 

NEWTONS 

4 

SPRING  RATE 

4.00000 

kN/m 

5 

END  COIL  NO. 

2.00000 

6 

SHEAR  FATIGUE  STRENGTH 

465.00000 

MPa 

7 

TENSILE  STRENGTH  CONSTANT 

1510.00000 

MPa 

8 

TENSILE  STRENGTH  WIRE  POWER 

-0.20100 

9 

SHEAR  MODULUS 

79000.00000 

MPa 

10 

WEIGHT  DENSITY 

76.50000 

kN/m**3 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 


LOW  HIGH  INITIAL 

1  WIRE  DIAMETER  1.00000  15.00000  5.00000  mm 

2  MEAN  COIL  DIAMETER  25.00000  500.00000  100.00000  mm 


THE 

7  CONSTRAINT  LIMITS  ARE: 

TYPE 

1 

ACTIVE  COIL 

0.00000 

NO. 

LOWER 

2 

FATIGUE  DESIGN  FACTOR 

1.50000 

LOWER 

3 

STATIC  DESIGN  FACTOR 

1.50000 

LOWER 

4 

SPRING  INDEX 

2.00000 

LOWER 

5 

MAXIMUM  SHEAR  STRESS 

0.00000 

MPa 

LOWER 

6 

OUTSIDE  COIL  DIAMETER 

0.00000 

mm 

LOWER 

7 

SPRING  WIRE  VOLUME 

0.00000 

mm**3 

LOWER 

MINIMIZE  THE  OBJECTIVE  FUNCTION, 

WHICH 

HAS  3  TERMS. 

OBJ  =  THE  LINEAR  SUM  OF: 

1 

WEIGHT 

IN 

NEWTONS 

TIMES 

0.0000 

2 

HEIGHT 

IN 

mm 

TIMES 

1.0000 

3 

SPRING  COIL  VOLUME 

IN 

mm**3 

TIMES 

0.0000 
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Table  15  Continued 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  53.5408  ,  ITS  COMPONENTS  ARE: 


1 

WEIGHT 

=  3.7045 

NEWTONS 

TIMES 

0.0000 

2 

HEIGHT 

=  53.541 

mm 

TIMES 

1.0000 

3 

SPRING  COIL 

VOLUME 

=0.17205E+07 

mm**3 

TIMES 

0.0000 

THE 

LAST  CHANGE 

IN  THE  OBJECTIVE 

FUNCTION 

=  0. 

381470E-04 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.625000E-02 


THE  7  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1  ACTIVE  COIL 

=  .53093 

NO. 

.00000 

LOWER 

2  FATIGUE  DESIGN  FACTOR 

=  2.2086 

1.5000 

LOWER 

3  STATIC  DESIGN  FACTOR 

=  1.5438 

1.5000 

LOWER 

4  SPRING  INDEX 

=  30.786 

2.0000 

LOWER 

5  MAXIMUM  SHEAR  STRESS 

=  312.12 

MPa 

.00000 

LOWER 

6  OUTSIDE  COIL  DIAMETER 

=  199.46 

mm 

.00000 

LOWER 

7  SPRING  WIRE  VOLUME 

=  47506. 

mm**3 

.00000 

LOWER 

Table  15  Continued 

Soring  Design  Output  File  For  Minimum  Height 


DESIGN  CHECK 


X{I) 


1  WIRE  DIAMETER  6.50000  mm 

2  MEAN  COIL  DIAMETER  220.00000  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  53.3470  ,  ITS  COMPONENTS  ARE: 


1  WEIGHT 

=  4.2351 

NEWTONS 

TIMES 

0.0000 

2  HEIGHT 

=  53.347 

mm 

TIMES 

1.0000 

3  SPRING  COIL  VOLUME 

=0.21495E+07 

mm**3 

TIMES 

0.0000 

THE  7  CONSTRAINT 

VALUES  ARE: 

LIMIT 

TYPE 

I  ACTIVE  COIL 

=  .41387  NO. 

.00000 

LOWER 

2  FATIGUE  DESIGN  FACTOR 

=  2.1515 

1.5000 

LOWER 

3  STATIC  DESIGN  FACTOR 

=  1.5020 

1.5000 

LOWER 

4  SPRING  INDEX 

=  33.846 

2.0000 

LOWER 

5  MAXIMUM  SHEAR  STRESS 

=  318.54  MPa 

.00000 

LOWER 

6  OUTSIDE  COIL  DIAMETER 

=  226.50  mm 

.00000 

LOWER 

7  SPRING  WIRE  VOLUME 

=  55361.  mm**3 

.00000 

LOWER 
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different  from  the  first.  The  spring  has  a  wire  diameter  of  6.3  mm  and  a  mean 
coil  diameter  of  196  mm  with  only  0.53  active  coils  and  a  spring  index  of  31. 
The  height  is  shorter  at  53.5  mm  but  the  weight  is  higher  at  3.7  Newtons  and 
the  spring  coil  volume  is  much  higher  at  1.72*10^  mm^.  This  design  was  also 
limited  by  the  static  design  factor. 

However,  it  was  reached  in  only  14  trials  and  the  static  design  factor 
was  slightly  higher  than  1.5  at  1.54.  These  two  facts  indicate  that  this  may 
not  be  an  absolute  minimum  height  design.  Changing  the  starting  position  as 
was  done  for  the  bushing  design  or  changing  the  sensitivity  ranges  on  the  two 
independent  design  parameters  would  give  the  optimizer  a  chance  to  search 
longer  and  find  a  slightly  better  optimal  design.  However,  since  we  were 
going  to  change  the  wire  size  and  mean  coil  diameter  to  standard  values 
anyway,  the  starting  position  and  sensitivity  ranges  were  not  changed.  The 
design  check  chosen  has  a  wire  diameter  of  6.5  mm  and  a  mean  coil  diameter  of 
220  mm  for  a  height  of  53.3  mm.  This  design  is  also  heavy  and  large  with  a 
fraction  of  an  active  coil  at  0.41  and  a  spring  index  of  34. 

Table  16  shows  the  output  file  for  a  design  which  minimizes  the  spring 
coil  volume.  The  only  changes  in  the  input  file  for  this  case  were  the 
weighting  coefficients  in  the  objective  function  list  to  select  spring  volume 
as  the  target  for  minimization  and  the  problem  title.  This  design  was 
achieved  in  28  steps  and  has  a  wire  diameter  of  2.2  mm  and  a  mean  coil 
diameter  of  7.5  mm.  The  spring  coil  volume  is  the  smallest  of  the  found 
designs  at  26,300  mm^  which  is  about  one-half  that  of  the  first  design. 
However  it  is  heavier  at  0.99  Newtons  and  much  longer  at  355  mm.  The  spring 
index  is  small  at  3.4  and  the  number  of  active  coils  is  large  at  140. 

Changing  the  wire  diameter  to  2.5  mm  and  the  mean  coil  diameter  to  11.5  mm 
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Table  16 


Design  Out.nut  File  For  Minimum  Coil  Volume 


HELICAL  COIL  SPRING  -  MINIMUM  COIL  VOLUME 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 

FIXED  DESIGN  REQUIREMENTS: 

1  MINIMUM  FORCE 

2  MAXIMUM  FORCE 

3  MINIMUM  FORCE  WHEN  SOLID 

4  SPRING  RATE 

5  END  COIL  NO. 

6  SHEAR  FATIGUE  STRENGTH 

7  TENSILE  STRENGTH  CONSTANT 

8  TENSILE  STRENGTH  WIRE  POWER 

9  SHEAR  MODULUS 
10  WEIGHT  DENSITY 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 

LOW  HIGH  INITIAL 


24.00000  NEWTONS 
120.00000  NEWTONS 
150.00000  NEWTONS 
4.00000  kN/m 
2.00000 

465.00000  MPa 
1510.00000  MPa 
-0.20100 

79000.00000  MPa 
76.50000  kN/ni**3 


1  WIRE  DIAMETER 

1.00000  15.00000  5.00000 

mm 

2  MEAN  COIL  DIAMETER 

25.00000  500.00000  100.00000 

mm 

THE 

7  CONSTRAINT  LIMITS  ARE: 

TYPE 

1 

ACTIVE  COIL 

0.00000 

NO. 

LOWER 

2 

FATIGUE  DESIGN  FACTOR 

1.50000 

LOWER 

3 

STATIC  DESIGN  FACTOR 

1 . 50000 

LOWER 

4 

SPRING  INDEX 

2.00000 

LOWER 

5 

MAXIMUM  SHEAR  STRESS 

0.00000 

MPa 

LOWER 

6 

OUTSIDE  COIL  DIAMETER 

0.00000 

mm 

LOWER 

7 

SPRING  WIRE  VOLUME 

0.00000 

mm**3 

LOWER 

MINIMIZE  THE  OBJECTIVE  FUNCTION,  WHICH  HAS  3  TERMS. 


OBJ  =  THE  LINEAR  SUM  OF: 

1  WEIGHT 

2  HEIGHT 

3  SPRING  COIL  VOLUME 


IN  NEWTONS  TIMES  0.0000 

IN  iTiin  TIMES  0.0000 

IN  niin**3  TIMES  1.0000 
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Table  16  Continued 


Spring  Design  Output  File  For  Minimum  Coil  Volume 


OPTIMIZATION  SUCCESSFUL  IN  28  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 


X(I) 

1  WIRE  DIAMETER 

2  MEAN  COIL  DIAMETER 


2.21384  mm 
7.51284  mm 


THE 

MINIMUM  OBJECTIVE  FUNCTION  = 

26353.6 

,  ITS  COMPONENTS  ARE: 

1 

WEIGHT 

=0.98587 

NEWTONS  TIMES 

0.0000 

2 

HEIGHT 

=  354.67 

mm  TIMES 

0.0000 

3 

SPRING  COIL  VOLUME 

=  26354. 

mm**3  TIMES 

1.0000 

THE 

LAST  CHANGE  IN  THE  OBJECTIVE 

FUNCTION 

=  -0.78I250E-02 

THE 

LAST  STEP  CHANGE  SIZE  FOR  THE 

DESIGN  VARIABLE  =  0.I95313E-03 

THE  7  CONSTRAINT 

VALUES  ARE: 

LIMIT 

TYPE 

1 

ACTIVE  COIL 

=  140.58 

NO. 

.00000 

LOWER 

2 

FATIGUE  DESIGN  FACTOR 

=  2.2583 

1.5000 

LOWER 

3 

STATIC  DESIGN  FACTOR 

=  1.5032 

1.5000 

LOWER 

4 

SPRING  INDEX 

=  3.3870 

2.0000 

LOWER 

5 

MAXIMUM  SHEAR  STRESS 

=  395.29 

MPa 

.00000 

LOWER 

6 

OUTSIDE  COIL  DIAMETER 

=  9.7067 

mm 

.00000 

LOWER 

7 

SPRING  WIRE  VOLUME 

=  12907. 

mm**3 

.00000 

LOWER 
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Table  16  Continued 

Soring  Design  Output  File  For  Minimum  Coil  Volume 


DESIGN  CHECK 

X(I) 

1  WIRE  DIAMETER  2.50000  mm 

2  MEAN  COIL  DIAMETER  11.50000  mm 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  31196.3 


ITS  COMPONENTS  ARE: 


1  WEIGHT 

2  HEIGHT 

3  SPRING  COIL  VOLUME 


=0.88738  NEWTONS  TIMES  0.0000 
=  202.66  mm  TIMES  0.0000 
=  31196.  mm**3  TIMES  1.0000 


THE  7  CONSTRAINT  VALUES  ARE: 


1 

ACTIVE  COIL 

=  63.408 

2 

FATIGUE  DESIGN  FACTOR 

=  2.2430 

3 

STATIC  DESIGN  FACTOR 

=  1.5367 

4 

SPRING  INDEX 

=  4.6000 

5 

MAXIMUM  SHEAR  STRESS 

=  377.29 

6 

OUTSIDE  COIL  DIAMETER 

=  14.000 

7 

SPRING  WIRE  VOLUME 

=  11600. 

LIMIT 

TYPE 

NO. 

.00000 

LOWER 

1.5000 

LOWER 

1.5000 

LOWER 

2.0000 

LOWER 

MPa 

.00000 

LOWER 

mm 

.00000 

LOWER 

mm**3 

.00000 

LOWER 
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with  the  design  check  provision  gives  a  design  with  a  spring  coil  volume  of 
31,200  mm^,  a  weight  of  only  0.89  Newtons  and  a  height  of  203  mm.  This  spring 
has  63.4  active  coils  with  a  spring  index  of  4.6. 

Figure  11  shows  all  three  optimum  designs  to  the  same  scale  for 
comparison.  In  Figure  10,  which  was  used  earlier  to  describe  the  design 
space,  three  crosses  mark  the  locations  of  the  three  optimal  designs  in  this 
design  space.  Each  design  is  at  or  near  the  static  strength  limit,  and  the 
minimum  weight  design,  with  its  cross  between  the  other  two,  is  the  design  we 
expect  from  practice.  It  is  the  most  economical  in  that  it  uses  the  least 
material  and  it  is  not  extremely  long  or  large  in  diameter.  However,  the 
other  designs  may  still  have  their  applications. 
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MINIMUM  HEIGHT 

OPTIMUM  SPRING  DESIGNS 
FIGURE  1  1 
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SPUR  GEAR  OPTIMIZATION 

The  third  example  problem  is  that  of  a  compact  gear  mesh,  which  is  to 
have  a  compact  size  at  a  given  life,  load,  reduction  and  speed.  For  steel 
involute  teeth,  the  loading  which  may  cause  failure  is  dynamic  due  to  the 
variation  in  contact  geometry  and  load  sharing  as  the  teeth  enter  and  leave 
the  gear  mesh.  Three  primary  modes  of  failure  are  possible:  1)  bending  tooth 
fracture,  2)  tooth  surface  pitting,  and  3)  tooth  tip  scoring.  When  fracture 
of  the  gear  teeth  due  to  bending  is  the  primary  mode  of  failure,  the  minimum 
number  of  teeth  which  avoids  interference  offers  the  strongest  gear  set  for  a 
given  size  [16].  However,  as  speeds  increase,  so  do  the  prospects  for  pitting 
and  scoring  modes  of  failure. 

For  a  given  combination  of  gear  material  and  lubrication  conditions,  the 
design  problem  can  be  formulated  in  terms  of  three  independent  design 
variables:  The  number  of  pinion  teeth,  Np;  the  diametral  pitch,  P^;  and  the 
gear  and  pinion  face  width,  f.  Although  many  designs  can  transmit  the  same 
power  at  the  same  input  and  output  speeds,  two  designs  will  be  sought.  The 
first  will  have  the  minimum  center  distance  between  the  input  and  output 
shafts  for  a  specified  reliability  life.  This  design  will  also  have  the 
minimum  gear  weight  and  volume.  The  second  will  have  the  maximum  life  for  a 
specified  center  distance. 

Theory 

All  acceptable  gear  designs  avoid  involute  interference.  For  equal 
pinion  and  gear  addenda,  involute  interference  will  occur  at  the  base  of  the 
pinion  tooth.  As  shown  in  Figure  12,  involute  interference  can  be  measured  by 
the  roll  angle,  6-^,  which  corresponds  to  the  distance  along  the  line  of  action 
from  the  base  circle  of  the  pinion  to  the  addendum  circle  of  the  gear. 


71 


INVOLUTE  INTERFERENCE  CHECK 
FIGURE  12 
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This  angle  must  be  positive  for  the  gear  tooth  tip  to  contact  the  pinion  tooth 
on  its  involute  surface  and  avoid  interference. 

In  performing  the  gear  tooth  load  analysis,  the  tangential  load  on  the 
gear  mesh,  F^,  is  the  pinion  torque  divided  by  the  pitch  radius  of  the  pinion. 
This  is  the  nominal  force  acting  between  the  gears.  The  force  along  the  line 
of  action  is  this  force  divided  by  the  cosine  of  the  pressure  angle,  0.  The 
pitch  line  velocity  is  the  rotational  speed  of  the  pinion  times  its  pitch 
radius,  which  is  the  speed  of  the  pitch  circles. 

To  estimate  the  dynamic  load,  one  can  use  the  AGMA  velocity  factor 
model  [17].  In  terms  of  a  gear  quality  number,  Q^,  the  AGMA  estimate  of  the 
sum  of  the  transmitted  load  and  the  dynamic  load  is: 


where 


A  =  50  + 


(28) 


In  equation  (28),  the  gear  quality  number,  Q^,  may  have  a  value  between  6  and 
11  with  11  corresponding  to  the  higher  quality  gear.  In  this  example,  all 
gear  stresses  and  lives  are  calculated  using  this  total  dynamic  load,  F^,  with 
a  quality  number,  =  9. 

As  noted  above,  gear  tooth  bending  fatigue,  gear  tooth  pitting  and  gear 
tip  scoring  are  the  three  most  probable  modes  of  failure  for  gear  teeth.  The 
bending  fatigue  model  uses  the  AGMA  J  factor  [17]  to  estimate  the  bending 
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stress  with  the  dynamic  load  at  the  highest  point  of  single  tooth  contact  on 
the  pinion,  which  produces  the  maximum  tooth  bending  stress.  The  formula  for 
the  bending  stress  is: 


F. 


(29) 


Figure  13  shows  the  dynamic  load  on  the  gear  tooth  and  the  Lewis  parabola 
which  describes  the  strongest  constant  strength  beam  inscribed  in  the  tooth. 

Gear  tooth  pitting  is  a  result  of  contact  stress  on  the  gear  tooth 
surface.  The  maximum  contact  stress  and  gear  tip  Hertzian  pressure  are 
calculated  [18]  as: 


n  f  cos  0 


1  /  Pi 


1  /  Pr 


1/2 


1  -  V/, 


1  -  V. 


(30) 


This  maximum  contact  stress  occurs  at  the  lowest  point  of  full  load  contact  on 
the  pinion  tooth.  Figure  14  shows  the  two  teeth  in  contact  at  this  point  and 
their  radii  of  curvature  on  which  this  stress  is  based.  The  small  parabola  on 
the  tooth  surface  shows  this  contact  pressure  distribution. 

The  gear  tip  Hertzian  pressure  uses  one-third  of  the  total  dynamic  load 
since  the  load  is  shared  unequally  between  two  tooth  pairs  at  this  point  due 
to  the  elastic  interaction  of  the  two  tooth  pairs  in  contact.  The  gear  tip 
scoring  model  includes  the  pressure  times  velocity  factor  and  the  critical  oil 
scoring  temperature  model  from  lubrication  theory.  The  normal  pressure  times 
sliding  velocity  is  proportional  to  the  frictional  power  loss  of  the  gear  set. 
This  factor  is  the  highest  for  contact  at  the  gear  tip,  with  the  normal 


74 


GEAR  TOOTH  BENDING  STRENGTH  MODEL 
FIGURE  13 
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PINION  AND  GEAR 
TOOTH  SURFACE  CURVATURE 


FIGURE  14 


pressure  equal  to  the  gear  tip  Hertzian  pressure.  The  sliding  velocity  at  the 
gear  tip  is  given  by: 


^s  "  ^2  *^a2  ■*'  ^a2^  ”  ^1 

The  lubricating  oil  flash  temperature  is  another  factor  used  to  monitor 
gear  tooth  scoring.  One  estimate  of  this  temperature  [18]  is  given  by: 

3/4 

I  Jfr  F.  1  I  0A5uX^.X^VV\ 

Tf  ^  Tg  A  1  I  -  1  (32) 


0.45//„/h^6^  1 

(Rj  + 


The  gear  life  model  is  based  on  surface  pitting  and  is  similar  to  the 
model  for  rolling  element  bearings  [19].  Its  reliability  estimation  is  based 
on  the  two-parameter  Weibull  distribution: 


1  1  1 

111  1^1 

Ln 

=  Ln  -  •  - 

(33) 

'  R  ' 

'  0.9  '  fjo 

reliability 

relationship  of  equation  (33) 

is  for  a  specific 

load  which  determines  the  fjQ  life.  This  load,  F,  is  related  to  the  component 
dynamic  capacity,  C,  as: 

*10  ■  I  ^  I 

where  the  dynamic  capacity  of  the  component,  C,  is  the  load  which  has  a 
90  percent  reliability  life  of  one  million  cycles. 

For  a  spur  gear  tooth,  the  dynamic  capacity,  C^,  can  be  expressed  as  a 
function  of  Buckingham’s  load-stress  factor  [20],  B,  which  is  a  material 
strength  constant: 
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r  =  fl  I  - ^ -  I  (35) 

^  '  l/pj  +  I/P2 

With  the  dynamic  capacity  expressed  in  this  form,  the  material  strength 
constant  serves  the  role  of  the  surface  fatigue  strength,  of  the  AGMA 

design  code.  A  relation  for  the  material  strength  constant  in  terms  of  the 
surface  fatigue  strength  is: 


The  dynamic  capacity  of  the  whole  gear  is  lower  than  that  of  a  single 
tooth  due  to  the  number  of  teeth  subjected  to  the  same  load: 


The  gear  and  pinion  weights  are  modeled  with  solid  discs  which  have 
radii  equal  to  the  pitch  radii  of  the  two  gears  and  thicknesses  equal  to  the 
face  width. 

Programming 

Table  17  summarizes  this  design  problem  in  terms  of  the  constants, 
design  parameters,  inequality  constraints  and  objective  functions.  There  are 
thirteen  constants,  three  design  parameters,  fourteen  constraints  and  two 
objective  functions  in  this  list. 

Included  in  the  list  of  constants  are:  Poisson’s  ratio  and  the  elastic 


modulus  for  the  gear  material:  the  nominal  pressure  angle,  0,  and  gear  ratio, 
n,  for  the  mesh;  the  transmitted  power,  Hp,  and  pinion  speed,  the  material 
weight  density,  y,  surface  constant,  B,  Weibull  slope,  b,  and  load-life 
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Table  17 


Sour  Gear  Mpsh  Optimization  Parameters 


Constants 

Design  Parameters 

Inequality  Constraints 

Objective  Function 

V 

Nl 

(Jj  <  0.001 

c  . 

mm 

E 

I’d 

f/Dj  <  0.5 

or 

0 

f 

C  >  0.0 

^^m^max 

n 

Wj  +  W2  >  0.0 

Hp 

Tj  >  0.0 

Wi 

>  0.0 

Y 

V  >  0.0 

B 

Fj  >  0.0 

b 

<  40 , 000 . 

P 

<  150,000. 

R 

a^t  <  150,000. 

h 

?  V  <  100. 

s 

7^  <  275. 

1  >  2.0 
m 
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exponent,  p;  the  desired  reliability  of  the  mesh,  R;  the  lubricant  base 
temperature,  T^;  and  the  surface  finish,  fj^. 

The  three  design  parameters  are:  1)  the  number  of  teeth  on  the  pinion, 
Nj',  the  diametral  pitch,  P^;  and  the  face  width,  f. 

Design  constraints  listed  are:  active  constraints  for  involute 
interference  and  the  width  to  diameter  ratio;  inactive  constraints  of  being 
positive  for  the  center  distance,  weight  sum,  input  torque,  transmitted  force, 
pitch  line  velocity,  and  dynamic  force;  active  bending  and  pitting  strength 
limits;  active  scoring  limits  of  pressure  times  velocity  and  flash 
temperature;  and  the  mesh  life  bound. 

The  objective  functions  include  the  center  distance  which  is  to  be 
minimized  and  the  mesh  life  which  is  to  be  maximized.  In  running  the  program, 
these  two  functions  will  not  be  active  in  the  same  case.  The  listed  values 
are  for  the  first  objective  function  of  minimum  center  distance.  When  the 
program  is  run  to  maximize  the  mesh  life,  the  center  distance  limit  is  changed 
to  an  active  limit  of  2.5  inches  and  the  mesh  life  constraint  is  changed  to  be 
positive  and  thus  inactive. 

Subroutines  BOUNDS  and  VALUES,  developed  for  this  example  gear  problem, 
are  listed  together  in  Table  18.  These  two  routines  call  a  series  of  analysis 
subroutines  to  evaluate  the  properties  of  the  gear  design.  These  subroutines: 
DYNAM,  MESH,  GLIFE,  GEARWT  and  TEMPER  are  also  listed  in  Table  18  along  with 
LEWIS  which  MESH  calls.  Subroutine  DYNAM  determines  the  dynamic  load  in  the 
mesh  with  the  AGMA  velocity  factor  calculation  using  equations  (27)  and  (28). 
Subroutine  MESH  performs  a  kinematic  analysis  of  the  gear  mesh  geometry  and 
calculates  the  bending  stresses  and  Hertzian  contact  stresses  on  the  teeth  in 
the  mesh  with  equations  (26)  and  (29)  through  (31).  MESH  also  calculates  the 
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Table  18 


Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 

SUBROUTINE  BOUNDS(CONST,NCO,X,NX,VCSTR,NCS) 

BOUNDS  DETERMINES  THE  PRESENT  CONSTRAINT 
FUNCTION  VALUES 

FOR  A  SPUR  GEAR  LIFE  DESIGN  EXAMPLE 


PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NCS  -  NUMBER  OF  INEQUALITY  CONSTRAINTS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  VARIABLES 

VCSTR  -  PRESENT  CONSTRAINT  VALUES 

X  -  PRESENT  DESIGN  VARIABLE  VALUES 


ALL  VALUES  ARE  IN  PROBLEM  UNITS  (UNSCALED) 


CONST(I)  =  mu 
CONST(2)  =  E 
CONST(3)  =  PHI 
CONST(4)  =  N 
CONST(5)  =  Hp 
CONST(6)  =  Np 
CONST(7)  =  GAMMA 
C0NST(8)  =  Be 
CONST(9)  =  b 
CONST{IO)=  p 
CONST(ll)=  Rel 
CONST(I2)=  Tb 
CONST ( 13 )=  RMS 


POISSON’S  RATIO 
ELASTIC  MODULUS  (PSI) 

PRESSURE  ANGLE  (DEGREES) 

GEAR  RATIO 

PINION  POWER  (HORSEPOWER) 

PINION  SPEED  (RPM) 

MATERIAL  WEIGHT  DENSITY  (LBS/IN**3) 

MATERIAL  STRENGTH  CONSTANT  (PSI) 

WEIBULL  SLOPE 

LOAD- LIFE  FACTOR 

DESIGN  RELIABILITY 

OIL  INLET  TEMPERATURE  (DEGREES  F) 

TOOTH  SURFACE  FINISH  (RMS) 


VCSTR(l)  = 
VCSTR(2)  = 
VCSTR(3)  = 
VCSTR(4)  = 
VCSTR(5)  = 
VCSTR(6)  = 
VCSTR(7)  = 
VCSTR(8)  = 
VCSTR(9)  = 
VCSTR(I0)= 
VCSTR(II)= 
VCSTR(12)= 

VCSTR(I3)= 

VCSTR(14)= 


DELI  -  INVOLUTE  INTERFERENCE  ANGLE  (RADIANS) 
LAMMIN  -  FACE  WIDTH  TO  PINION  DIAMETER  RATIO 
C  -  CENTER  DISTANCE  (INCHES) 

WT  -  PINION  &  GEAR  WEIGHTS  (POUNDS) 

Tq  -  PINION  TORQUE  (POUND  -  INCHES) 

Ft  -  TRANSMITTED  FORCE  (POUNDS) 

V  -  PITCH  LINE  VELOCITY  (FT/MIN) 

Fd  -  DYNAMIC  LOAD  (POUNDS) 

BDSTR  -  AGMA  BENDING  STRESS  (PSI) 

HZSTR  -  CONTACT  PRESSURE  (PSI) 

TIPPR  -  GEAR  TIP  HERTZ  CONTACT  PRESSURE  (PSI) 
PVF  -  PRESSURE  TIMES  VELOCITY  FACTOR 

(10**6  PSI-FT/MIN) 

Tf  -  FLASH  TEMPERATURE  (DEGREES  F) 

Lm  -  MESH  LIFE  (10**3  HOURS) 


81 


o  o  o  o  o 


Table  18  Continued 


Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 


X(l)  =  nl 
X(2)  =  Pd 
X(3)  =  Wd 


-  NUMBER  OF  PINION  TEETH 

-  DIAMETRAL  PITCH  (1.0/IN) 

-  FACE  WIDTH  (INCHES) 


DIMENSION  X(NX) ,VCSTR(NCS) ,CONST(NCO) 

REAL  LM,LP,LG 
PI  =  3.14159265 
POI  =  CONST(l) 

E  =  CONST(2) 

PHI  =  CONST(3) 

RATIO  =  CONST(4) 

HP  =  CONST(5) 

RPM  =  CONST(6) 

GAM  =  CONST(7) 

BC  =  C0NST(8) 

B  =  CONST(9) 

P  =  CONST(IO) 

REL  =  CONST(ll) 

TB  =  CONST(12) 

RMS  =  CONST(13) 

TNI  =  X(l) 

PD  =  X(2) 

WD  =  X(3) 

TN2  =  RATIO*TNl 
PHIR  =  PHI*PI/180.0 
RP  =  0.5*TN1/PD 
RG  =  RATIO*RP 
C  =  RP*(  1.0  +  RATIO  ) 

V  =  RPM*RP*PI/6.0 
TQ  =  HP*63025.0/RPM 
FT  =  TQ/RP 
CALL  DYNAM(FT,V,DL) 

TDQ  =  RP*DL 

CALL  MESH ( PH  I , PD , TN 1 , TN2 , E , PO I , WD , TDQ , RPM , CMP , DEL 1 , RH 1 , 
1  PVF,TIPHZJIPBS,HZSTR,BDSTR) 

RH2  =  C*SIN(PHIR)  -  RHl 

CALL  GLIFE(DL,WD,RH1,RH2,TN1,TN2,RPM,BC,B,P,REL, 

1  I  M  I  P  I G  CYLP  CYLGi 
C A  L  L ’ G  EARwf ( GAM ’ GAM , R  P , RG , WD , PW , G W ) 

CALL  TEMPER(PD,PHI,TN1,TN2,WD,DL,V,TB,RMS,TF,CSN) 
VCSTR(l)  =  DELI 
VCSTR(2)  =  WD/(2.0*RP) 

VCSTR(3)  =  C 
VCSTR(4)  =  PW  +  GW 
VCSTR(5)  =  TQ 
VCSTR(6)  =  FT 
VCSTR(7)  =  V 
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Table  18  Continued 

Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 


VCSTR(8)  =  DL 
VCSTR(9)  =  BDSTR 
VCSTR(IO)  =  HZSTR 
VCSTR(ll)  =  TIPHZ 
VCSTR(12)  =  PVF/1000000.0 
VCSTR(13)  =  TF 
VCSTR(14)  =  LM/1000.0 
RETURN 
END 
C 

SUBROUTINE  VALUES(CONST,NCO,X, NX, OBJECT, NOB) 

VALUES  DETERMINES  THE  PRESENT  DESIGN 
OBJECTIVE  FUNCTION  VALUES 


PARAMETERS: 

CONST  -  FIXED  DESIGN  CONSTANT 

NCO  -  NUMBER  OF  DESIGN  CONSTANTS 

NOB  -  NUMBER  OF  OBJECTIVE  FUNCTION  TERMS 

NX  -  NUMBER  OF  INDEPENDENT  DESIGN  VARIABLES 

OBJECT  -  PRESENT  OBJECTIVE  FUNCTION  VALUES 
X  -  PRESENT  DESIGN  VARIABLE  VALUES 

ALL  VALUES  ARE  IN  PROBLEM  UNITS  (UNSCALED) 


CONST(l)  =  mu 
CONST{2)  =  E 
CONST(3)  =  PHI 
CONST(4)  =  N 
CONST(5)  =  Hp 
CONST(6)  =  Np 
CONST (7)  =  GAMMA 
C0NST(8)  =  Be 
CONST(9)  =  b 
CONST(10)=  p 
CONST(ll)=  Rel 
CONST(12)=  Tb 
CONST(13)=  RMS 


-  POISSON’S  RATIO 

-  ELASTIC  MODULUS  (PSI) 

-  PRESSURE  ANGLE  (DEGREES) 

-  GEAR  RATIO 

-  PINION  POWER  (HORSEPOWER) 

-  PINION  SPEED  (RPM) 

-  MATERIAL  WEIGHT  DENSITY  (LBS/IN**3) 

-  MATERIAL  STRENGTH  CONSTANT  (PSI) 

-  WEIBULL  SLOPE 

-  LOAD- LIFE  FACTOR 

-  DESIGN  RELIABILITY 

-  BASE  TEMPERATURE  (DEG  FAHR) 

-  TOOTH  SURFACE  FINISH  (RMS) 


OBJECT(l)  =  C  -  CENTER  DISTANCE  (INCHES) 

OBJECT(2)  =  Lm  -  MESH  LIFE  (10**3  HOURS) 


X(l)  =  nl 
X(2)  =  Pd 
X(3)  =  Wd 


NUMBER  OF  PINION  TEETH 
DIAMETRAL  PITCH  (1.0/IN) 
FACE  WIDTH  (INCHES) 
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Table  18  Continued 


Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 


DIMENSION  X(NX) ,0BJECT(N0B) ,C0NST(NC0) 

REAL  LM,LP,LG 
PI  =  3.14159265 
POI  =  CONST(l) 

E  =  C0NST(2) 

PHI  =  C0NST(3) 

RATIO  =  C0NST(4) 

HP  =  C0NST(5) 

RPM  =  C0NST(6) 

BC  =  C0NST(8) 

B  =  C0NST(9) 

P  =  CONST(IO) 

REL  =  CONST(ll) 

TNI  =  X(l) 

PD  =  X(2) 

WD  =  X(3) 

TN2  =  RATIO*TNl 
PHIR  =  PHI*PI/180.0 
RP  =  0.5*TN1/PD 
C  =  RP*(  1.0  +  RATIO  ) 

V  =  RPM*RP*PI/6.0 
TQ  =  HP*63025.0/RPM 
FT  =  TQ/RP 
CALL  DYNAM(FT,V,DL) 

TDQ  =  RP*DL 

CALL  MESH(PHI,PD,TN1,TN2,E,P0I,WD,TDQ,RPM,CMP,DEL1,RH1, 

1  PVFJIPHZ,TIPBS,HZSTR,BDSTR) 

RH2  =  C*SIN(PHIR)  -  RHl 

CALL  GLIFE(DL,WD,RH1,RH2,TN1JN2,RPM,BC,B,P,REL, 

1  LM,LP,LG,CYLP,CYLG) 

OBJECT(l)  =  C 
OBJECT(2)  =  LM/1000.0 
RETURN 
END 

SUBROUTINE  MESH(PHI,PD,TN1JN2,EL,P0I,WD,TQ, RPM, CMP, DELI, ROB, 
1  PVF,SQA,AGBN2,SQB,AGBN1) 

MESH 

THIS  ROUTINE  CALCULATES  THE  GEAR  TOOTH  CONTACT  GEOMETRY 
INPUTS: 

PHI  -  NOMINAL  PRESSURE  ANGLE  (DEGREES) 

PD  -  DIAMETRAL  PITCH  (1.0/INCH) 

TNI  -  NUMBER  OF  TEETH  ON  PINION 

TN2  -  NUMBER  OF  TEETH  ON  GEAR 
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Table  18  Continued 

Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 

EL  -  ELASTIC  MODULUS  OF  TEETH  (PSI) 

POI  -  POISSON’S  RATIO  FOR  TEETH 
WD  -  TOOTH  FACE  WIDTH  (INCHES) 

TQ  -  TORQUE  APPLIED  TO  PINION  (POUND  -  INCHES) 

RPM  -  SPEED  OF  PINION  IN  REVOLUTIONS  PER  MINUTE 

OUTPUTS: 

CMP  -  MESH  CONTACT  RATIO 

DELI  -  DISTANCE  FROM  PINION  BASE  CIRCLE  TANGENCY 

TO  TIP  OF  GEAR  TOOTH  ALONG  LINE  OF  ACTION 
MEASURED  AS  A  PINION  ROLL  ANGLE 
ROB  -  PINION  TOOTH  SURFACE  RADIUS  OF  CURVATURE 

AT  THE  LOWEST  POINT  OF  SINGLE  TOOTH  CONTACT  C 
PVF  -  SCORING  FACTOR  (PSI-FT/MIN  /  10**6) 

SQA  -  HERTZIAN  CONTACT  STRESS  AT  TIP  OF  GEAR  TOOTH 

(PSI) 

AGBN2  -  AGMA  BENDING  STRESS  FOR  HALF  LOAD  AT  TIP 
OF  GEAR  TOOTH  (PSI) 

SQB  -  HERTZIAN  CONTACT  STRESS  AT  LOWEST  POINT  OF 
SINGLE  TOOTH  CONTACT  (PSI) 

AGBNl  -  AGMA  BENDING  STRESS  FOR  FULL  LOAD  AT 

HIGHEST  POINT  OF  SINGLE  TOOTH  CONTACT  (PSI)  C 
ACS (A) =ATAN (SQRT( 1 . -A*A) /A) 

ASN ( A) =ATAN ( A/SQRT ( 1 . -A*A) ) 

FINV(A)=SIN(A)/COS(A)-A 

NW=0 

PI=3. 14159265 

IT  IS  ASSUMED  THAT  THE  PINION  IS  THE  DRIVING  GEAR 


J=I 
FJ=J 
ADD=1 . 

DED=1.25 

RFR=0.3 

IF(Pd!gE.20.0)  DED=1.2+0.002*PD 

ADD1=ADD 

ADD2=ADD 

DDD1=DED 

AD1=ADD1/PD 

AD2=ADD2/PD 

DD1=DDD1/PD 

RF=RFR/PD 

PC=PI/PD 

RAD=180./PI 

P=PHI/RAD 

COP=COS(P) 

SIP=SIN(P) 
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Table  18  Continued 

Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 


PB=PC*C0P 

R1=0.5*PC*TN1/PI 

R2=0.5*PC*TN2/PI 

RB1=R1*C0P 

RB2=R2*COP 

RA1=R1+AD1 

RA2=R2+AD2*FJ 

C=FJ*R1+R2 

C1=SQRT(RA1*RA1-RB1*RB1) 

C2=SQRT(RA2*RA2-RB2*RB2) 

Z=C1+(C2-C*SIP)*FJ 

CMP=Z/PB 

AA1=(C2-R2*SIP)/RB1*FJ 

AA2=(C2-R2*SIP)/RB2*FJ 

DEL1=SIP/C0P-AA1 

IF  DELI  LESS  THAN  0  INTERFERENCE  WILL  HAPPEN 

SET  EQUAL  TO  A  SMALL  VALUE  TO  GENERATE 
LARGE  CONTACT  STRESSES 

IF(DELl.LT.O.O)  DELI  =  0.001*SIP/COP 
AB2=(C1-R1*SIP)/RB2 
DEL2=SIP/COP-AB2 
IF(J.LT.O)  DEL2=SIP/COP-AA2 

IF  DEL2  LESS  THAN  0  INTERFERENCE  WILL  HAPPEN 

SET  EQUAL  TO  A  SMALL  VALUE  TO  GENERATE 
LARGE  CONTACT  STRESSES 

IF(DEL2.LT.0.0)  DEL2  =  0.001*SIP/COP 

BL1=(Z-PB)/RB1 

BH1=(2.*PB-Z)/RB1 

TL1=DEL1+BL1 

TU1=TL1+BH1 

RAT=TN2/TN1 

IF(J.GT.O)  GO  TO  5 

TP=FINV(P) 

PA1=ACS(RB1/RA1) 

PA2=ACS(RB2/RA2) 

TA1=FINV(PA1) 

TA2=FINV(PA2) 

BETA1=ACS ( ( RA2*RA2-RA1*RA1-C*C) / ( 2 . 0*C*RA1 ) ) 
BETA2=ASN(RA1*SIN(BETA1)/RA2) 

GAMMA1=BETA1+TA1-TP 

GLIM2=BETA2+TA2-TP 

GAMMA2=TN1*GAMMA1/TN2 
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Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 


IF(GAMMA2.GT.GLIM2)  GO  TO  5 
XINT=RA2*(GLIM2-GAMMA2) 

WRITE  (NW,122)  XINT 

122  F0RMAT(//’  SECONDARY  INVOLUTE  INTERFERENCE’/ 

1  ’  TOOTH  TIP  OVERLAP  =  ’,F8.5,’  INCHES  AT  ENTRY’/ 

2  ’  INCREASE  N2/N1  TO  AVOID  THIS  INTERFERENCE’/) 

5  CONTINUE 

FSIN  =  TQ/Rl 
FTIP  =  FSIN/3.0 
EO=EL/(I.-POI*POI) 

R0A=RB1*DELI 

R0B=RBI*TL1 

CAP1=1./(R0A*(1.-FJ*R0A/(C*SIP))) 

CAP2= 1 . / ( ROB* ( 1 . - F J*ROB/ ( C*S I P ) ) ) 

SBA=SQRT(2.*FTIP/(COP*PI*EO*WD*CAPl)) 

SQA=FTIP/(COP*PI*WD*SBA) 

SBB=SQRT(2 .*FSIN/ (C0P*PI*E0*WD*CAP2) ) 
SQB=FSIN/(COP*PI*WD*SBB) 

PU1=ATAN(TU1)*RAD 

TC=PC/2. 

CALL  LEWIS(PHI,DD1,TC,RF,PD,TN1,PUI,YS,AKS) 
SBND1=FSIN*PD/(WD*YS) 

AGBN1=SBNDI*AKS 

TT1=TU1+BL1 

PT1=ATAN(TT1)*RAD 

CALL  LEWIS(PHI,DD1,TC,RF,PD.TN1,PT1,YT,AKT) 
SBND2=FTIP*PD/(WD*YT) 

AGBN2=SBND2*AKT 

RTP1=RA2*SIN(AA2)/SIN(AA1) 

OMEGA=RPM*2.*PI 

VLG= (OMEGA/RAT) *RA2*S IN ( P+AA2 ) 

VLP=OMEGA*RTPl*S IN ( P-AAl ) 

PVF=SQA*(VLG-VLP)/12. 

RETURN 

END 

C 

SUBROUTINE  LEWIS(PHI,B,TC,RF,PDJN,PHIA,Y,AK) 

LEWIS  CALCULATES  THE  LEWIS  FORM  FACTOR  FOR  A  PINION  TOOTH 

PHI  =  THE  PITCH  LINE  PRESSURE  ANGLE  IN  DEGREES 
B  =  THE  PINION  DEDENDUM  IN  INCHES 
TC  =  THE  PITCH  CIRCLE  TOOTH  THICKNESS  IN  INCHES 
RF  =  THE  CUTTER  TIP  RADIUS  IN  INCHES 
PD  =  THE  DIAMETRAL  PITCH  IN  1. /INCHES 
TN  =  THE  PINION  TOOTH  NUMBER 
PHIA  =  THE  PRESSURE  ANGLE  AT  THE  POINT  OF  CONTACT 
FOR  WHICH  Y  IS  CALCULATED  IN  DEGREES 
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Y  =  THE  LEWIS  FORM  FACTOR 
AK  =  THE  TOOTH  FILLET  STRESS  CONCENTRATION  FACTOR 

FINV(X)=SIN(X)/COS(X)-X 

PI=3. 14159265 

NW  =  0 

ACC=.001 

P=PHI*PI/180.0 

PHA=PHIA*PI/180.0 

RP=TN/(2.*PD) 

YB=B-RF 

C0P=C0S{P) 

SIP=SIN{P) 

DELTA=PI/{2.*PD)-TC/2.0-YB*SIP/C0P-RF/C0P 

RC=RP*C0P/(C0S(PHA-TC/(2.0*RP)-FINV(P)+FINV(PHA))) 

BETA=PI/TN-DELTA/RP 

ER=1. 

DTH=.01 

THETA=.00 

DO  3  1=1,500 

THETA=THETA+DTH 

XB=RP*THETA 

RB=SQRT(YB**2+XB**2) 

COA=COS(BETA+THETA) 

SIA=SIN(BETA+THETA) 

XE=(RP-YB)*SIA-XB*COA-(RF/RB)*(YB*SIA+XB*COA) 
YE=(RP-YB)*COA+XB*SIA+(RF/RB)*(XB*SIA-YB*COA) 
SLOPE=-((l.+(YB/XB)*(SIA/COA))/(YB/XB-SIA/COA)) 
ERR=SLOPE+(2.*(RC-YE))/XE 
IF(ABS(ERR).LT.ACC)GO  TO  4 
I F ( ERR/ ER . LT . 0 . ) DTH=-DTH/2 . 

ER=ERR 

3  CONTINUE 
WRITE  (NW,1) 

1  F0RMAT(’0  ITERATION  FOR  THETA  UNSUCCESSFUL’) 

Y=0. 

AK=1.0 
GO  TO  5 

4  Y=2.*((XE**2)/(RC-YE))*PD/3. 

AP=0.4583662*P 

AH=0.340-AP 

AL=0.316-AP 

AM=0.290+AP 

RAP=RF+YB*YB/(RP+YB) 

TR=(2.0*XE/RAP)**AL 

TY=(2.0*XE/BH)**AM 

AK=AH+TR*TY 
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Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 

5  RETURN 
END 

SUBROUTINE  GLIFE(FORCE,WD,RHP,RHG,TNP,TNG,PINRPM,BC,B,P,REL, 

1  LM,LP,LG,CYLP,CYLG) 

SUBROUTINE  TO  CALCULATE  THE  LIFE  OF  A  GEAR  MESH 

INPUTS 

FORCE  =  NORMAL  TOOTH  LOAD  IN  POUNDS 

WD  =  EFFECTIVE  FACE  WIDTH  IN  INCHES 

RHP  =  PINION  CONTACT  RADIUS  OF  CURVATURE  IN  INCHES 

RHG  =  GEAR  CONTACT  RADIUS  OF  CURVATURE  IN  INCHES 

TNP  =  NUMBER  OF  TEETH  ON  PINION  AS  A  REAL  VARIABLE 

TNG  =  NUMBER  OF  TEETH  ON  GEAR  AS  A  REAL  VARIABLE 

PINRPM  =  PINION  SPEED  IN  REVOLUTIONS  PER  MINUTE 

BC  =  MATERIAL  STRENGTH  CONSTANT  IN  PS I 

B  =  WEIBULL  SLOPE 

P  =  LOAD  LIFE  FACTOR 

REL  =  DESIRED  RELIABILITY  AS  A  DECIMAL 

OUTPUTS 

LM  =  MESH  LIFE  IN  HOURS 

LP  =  PINION  LIFE  IN  HOURS 

LG  =  GEAR  LIFE  IN  HOURS 

CYLP  =  PINION  LIFE  IN  MILLION  ROTATIONS 

CYLG  =  GEAR  LIFE  IN  MILLION  ROTATIONS 


REAL  LM,LP,LG 

CT  =  BC*WD/( 1.0/RHP  +  1.0/RHG) 

EX  =  1.0/(B*P) 

CP  =  CT*(1.0/TNP)**EX 
CG  =  CT*(1.0/TNG)**EX 
PLIO  =  (CP/FORCE)**P 
GLIO  =  (CG/FORCE)**P 
BR  =  1.0/B 

RATLF  =  (AL0G(REL)/AL0G(0.9))**BR 

CYLP  =  PL10*RATLF 

CYLG  =  GL10*RATLF 

LP  =  CYLP*1000000. 0/(60. 0*PINRPM) 

LG  =  CYLG*1000000.0*TNG/(60.0*PINRPM*TNP) 
LM  =  1.0/((1.0/TP)**B  +  (1.0/TG)**B)**BR 
RETURN 
END 


89 


oooooooooooooooooo  o  ooooooooooooo 


Table  18  Continued 


Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 
SUBROUTINE  DYNAM  (FT,V,DL) 

DYLOAD  CALCULATES  THE  DYNAMIC  LOAD  OF  THE  MESHING  OF 
THE  GEARS  USING  THE  AGMA  VELOCITY  FACTOR 

DL  =  THE  DYNAMIC  LOAD  DUE  TO  THE  MESHING 

INPUTS: 

FT  =  THE  TRANSMITTED  LOAD  (POUNDS) 

V  =  THE  PITCH  LINE  VELOCITY  (FT/MIN) 

OUTPUT: 

DL  =  THE  TOTAL  DYNAMIC  LOAD  (POUNDS) 

QV  =  9.0 

A  =  50.0  +  56.0*(1.0  -  (12.0  -  QV)**(2.0/3.0))/4.0 

VF  =  (  A  +  SQRT(V)  )/A 

DL  =  FT*VF 

RETURN 

END 

SUBROUTINE  TEMPER(PD,PHI,TN1,TN2,WD,DL,V,TB,RMS,TF,CSN) 

TEMPER  IS  A  SUBROUTINE  THAT  CALCULATES  THE  FLASH 
TEMPERATURE,  THE  CRITICAL  SCORING  NUMBER  AND  THE  MINIMUM 
ELASTOHYDRODYNAMIC  FILM  THICKNESS  FOR  TWO  GEARS  IN  MESH 


PHI 

5= 

PRESSURE  ANGLE  IN  DEGREES 

TN1,TN2 

= 

NUMBER  OF  TEETH  ON  PINION  AND  GEAR 

PD 

= 

DIAMETRAL  PITCH 

V 

= 

PITCH  LINE  VELOCITY  (FEET/MINUTE) 

RMS 

= 

TOOTH  SURFACE  FINISH  (RMS) 

RA2 

= 

GEAR  ADDENDUM  RADIUS 

RB2 

GEAR  BASE  RADIUS 

WD 

FACE  WIDTH 

DL 

TOOTH  DYNAMIC  LOAD  (POUNDS) 

TB 

= 

INLET  OIL  TEMPERATURE 

C2 

= 

RADIUS  OF  CURVATURE  AT  TIP  OF  GEAR  TOOTH 

TF 

= 

FLASH  TEMPERATURE 

FT 

NORMAL  TOOTH  LOAD  IN  POUNDS 

PI  =  3.14159265 

P  =  PHI*PI/180.0 
COP  =  COS(P) 

SIP  =  SIN(P) 

R1  =  TN1/(2.0*PD) 
R2  =  TN2/(2.0*PD) 
C  =  R1  +  R2 
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ooooooooo  o  ooo  ooo  ooo 


Table  18  Continued 


Sour  Gear  Mesh  Constraint  and  Objective  Function  Evaluation  Routines 

GR  =  TN2/TN1 

RPM  =  6.0*V/(PI*R1*C0P) 

RAl  =  R1  +  1.0/PD 
RBI  =  R1*C0P 

Cl  =  SQRT(RA1*RA1  -  RB1*RB1) 

TAT  =  Cl/RBl 
TAP  =  SI P/COP 
GY  =  TAT/TAP  -  1.0 

FD  IS  THE  NORMAL  FORCE  ON  THE  GEAR 

FD  =  DL/3.0 

CRITICAL  SCORING  TEMPERATURE  OF  THE  GEAR 


XM  =  1.75 

AX  =  ABS(SQRT(1.0  +  GY)  -  SQRT(1.0  -  GY/GR)) 

BX  =  ((1.0  +  GY)*(GR  -  GY))**0.25 
XG  =  0.51*SQRT(GR  +  1.0)*AX/BX 
SW  =  (FD/WD)**(0.75) 

TI  =  SW*(0. 027*50.0/(50.0  -  RMS))*XM*XG*SQRT(V)/(C**0.25) 

TF  =  TB  +  TI 

CRITICAL  SCORING  NUMBER 

CSN  =  SW*SQRT(RPM)/(PD**.25) 

RETURN 

END 

SUBROUTINE  GEARWT  (RH01,RH02,R1,R2,WD,PW,GW) 

GEARWT  CALCULATES  THE  WEIGHT  OF  THE  PINION  AND  GEAR 
IN  MESH 

RH01,RH02  =  MATERIAL  DENSITY  OF  THE  RESPECTIVE  GEARS 

R1,R2  =  PITCH  RADIUS  OF  THE  PINION  AND  GEAR 

WD  =  FACE  WIDTH  OF  GEARS 

PW  =  PINION  WEIGHT 

GW  =  GEAR  WEIGHT 


PI  =  3.14159265 

PW  =  RH01*PI*R1**2  *WD 

GW  =  RH02*PI*R2**2  *WD 

RETURN 

END 
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pressure  times  velocity  at  the  gear  tip  to  model  potential  scoring.  In  order 
to  calculate  the  bending  stresses,  subroutine  MESH  calls  subroutine  LEWIS 
which  performs  an  interval  halving  iteration  to  determine  the  size  of  the 
largest  inscribed  parabola  in  the  tooth.  Subroutine  GLIFE  performs  a  Weibull 
analysis  with  equations  (33)  through  (37)  to  estimate  the  life  of  the  gear 
mesh  for  the  reliability  specified  in  the  constant  list.  Subroutine  GEARWT 
estimates  the  weights  of  the  pinion  and  the  gear  and  subroutine  TEMPER 
calculates  the  flash  temperature  for  the  mesh  using  equation  (32)  to  determine 
the  potential  for  scoring  damage. 

Any  permissible  subroutine  names  may  be  used  in  programming  for  a 
specific  problem  as  long  as  the  names  of  the  subroutines  in  SEEK  are  avoided. 
These  names  are  listed  in  the  file  SEEK.DOC,  which  is  Appendix  A,  and  in 
Table  1  of  the  programming  section  and  are:  BACK,  BOUNCE,  BOUNDS,  CHECK, 
GRADNT,  MERIT,  RESIZE,  SCAN,  SCOUT,  SIZE,  UNIT,  VALUES  and  WALL.  Once  the 
problem  subroutines  have  been  written,  they  need  to  be  compiled  and  linked  to 
the  compiled  program  SEEK. OBJ  to  produce  an  executable  program  with  a  name 
such  as  GEAR.EXE.  With  this  program,  various  optimal  gear  designs  can  be 
found  using  different  input  files. 

Numerical  Example 

For  an  example,  let  us  consider  the  design  of  a  gear  set  to  transmit 
10  horsepower  from  a  shaft  turning  at  4,500  RPM  to  an  output  shaft  turning  at 
3,000  RPM.  The  center  distance  of  the  gears  should  be  minimized  for  a  mesh 
life  of  at  least  2,000  hours  with  a  reliability  of  90  percent.  The  face  width 
ratio  is  to  be  less  than  or  equal  to  0.5,  the  material  strengths  are 
40,000  psi  in  bending  and  150,000  psi  surface  endurance,  the  PV  factor  limit 
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is  100  million  psi-ft/min  and  the  flash  temperature  limit  is  275  degrees 
Fahrenheit. 

Table  19  is  a  listing  of  the  input  file  defining  this  example.  Rules 
for  writing  the  input  data  file  are  given  in  "SEEK.DOC,"  which  is  listed  in 
Appendix  A.  The  data  file  starts  with  the  problem  title.  On  the  next  line  is 
the  number  thirteen,  and  the  following  thirty-nine  lines  contain  the  problem 
constants  as  listed  in  Table  17,  their  definitions  and  their  units.  The 
material  surface  constant  of  9,800  psi  corresponds  to  a  surface  compression 
endurance  strength  of  200,000  psi  at  10^  fatigue  cycles  and  a  reliability  of 
90  percent.  The  load-life  factor  of  8.93  is  from  the  ANSI/A6MA  2001  888 
Standard,  and  the  Weibull  slope  of  2.5  is  from  the  NASA  Lewis  gear  test  data 
[21].  Following  the  constants  is  a  line  with  the  number  three  for  the  number 
of  independent  design  parameters,  which  are  on  the  next  nine  lines  with  their 
ranges,  initial  values,  definitions  and  units.  The  number  fourteen  follows  on 
the  next  line,  with  the  fourteen  constraint  limits,  their  directions  and  their 
names  and  units  on  the  following  forty-two  lines.  After  these  comes  a  line 
with  the  optimization  direction,  MIN,  and  a  line  with  the  number  two  for  the 
number  of  objective  function  terms.  The  first  term  is  the  center  distance 
which  has  a  weighting  coefficient  of  one  and  the  units  of  inches,  and  the 
second  is  the  mesh  life  with  a  weighting  coefficient  of  zero  and  units  of 
thousand  hours. 

Running  the  program  GEAR.EXE  with  this  data  file  produced  the  output 
data  file  listed  in  Table  20.  The  first  page  lists  the  problem  constants, 
design  parameters  and  constraint  limits  as  provided  by  the  input  data  file  but 
in  a  little  more  readable  form.  The  second  page  lists  the  merit  function 
terms  and  notes  that  an  optimum  was  found  in  thirty-three  steps.  It  then 


93 


Table  19 


Sour  Gear  Mesh  Input  File  for  Minimum  Size 

COMPACT  SPUR  GEAR  WITH  A  REQUIRED  LIFE 
13 

0.25 

POISSON’S  RATIO 

30000000.0 
ELASTIC  MODULUS 
PS  I 
20.0 

PRESSURE  ANGLE 
DEGREES 

1.5 

GEAR  RATIO 

10.0 

TRANSMITTED  POWER 

HORSEPOWER 

4500.0 

PINION  SPEED 

RPM 

0.283 

MATERIAL  WEIGHT  DENSITY 

LBS/IN**3 

9800.0 

MATERIAL  SURFACE  CONSTANT 
PSI 

2.5 

WEIBULL  SLOPE 
8  93 

LOAD-LIFE  FACTOR 
0.90 

RELIABILITY 

120.0 

BASE  TEMPERATURE 
DEGREES  F 
32 

TOOTH  SURFACE  FINISH 
RMS 
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Table  19  Continued 


Sour  Gear  Mesh  Input  File  for  Minimum  Size 


3 

10.0  100.0  40.0 
PINION  TEETH 

4.0  28.0  10.0 
DIAMETRAL  PITCH 
1.0/INCHES 
0.5  5.0  2.5 
FACE  WIDTH 
INCHES 
14 

LOWER  0.001 
INVOL.  INTERFERENCE 
RADIANS 
UPPER  0.5 

FACE  WIDTH  TO  DIAMETER 

RATIO 

LOWER  0.0 

CENTER  DISTANCE 

INCHES 

LOWER  0.0 

GEAR  AND  PINION  WEIGHT 

POUNDS 

LOWER  0.0 

PINION  TORQUE 

LB-IN 

LOWER  0.0 

TRANSMITTED  FORCE 

POUNDS 

LOWER  0.0 

PITCH  LINE  VELOCITY 

FT/MIN 

LOWER  0.0 

TOTAL  DYNAMIC  LOAD 
POUNDS 

UPPER  40000.0 
AGMA  BENDING  STRESS 
PS  I 

UPPER  150000.0 

FULL  LOAD  CONTACT  STRESS 

PS  I 

UPPER  150000.0 

GEAR  TIP  HERTZ  PRESSURE 

PSI 

UPPER  100.0 
PV  FACTOR 
M  PSI-FT/MIN 
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Table  19  Continued 


Sour  Gear  Mesh  Input  File  for  Minimuin  Size 


UPPER  275.0 
FLASH  TEMPERATURE 
DEGREES  F 
LOWER  2.0 
MESH  LIFE 
10**3  HOURS 
MIN 
2 

1.0 

CENTER  DISTANCE 
INCHES 
0.0 

MESH  LIFE 
10**3  HOURS 
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Table  20 


Sour  Gear  Mesh  Output  File  for  Minimum  Size 
COMPACT  SPUR  GEAR  WITH  A  REQUIRED  LIFE 


DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 


FIXED  DESIGN  REQUIREMENTS: 


1  POISSON’S  RATIO 

2  ELASTIC  MODULUS 

3  PRESSURE  ANGLE 

4  GEAR  RATIO 

5  TRANSMITTED  POWER 

6  PINION  SPEED 

7  MATERIAL  WEIGHT  DENSITY 

8  MATERIAL  SURFACE  CONSTANT 

9  WEIBULL  SLOPE 

10  LOAD-LIFE  FACTOR 

11  RELIABILITY 

12  BASE  TEMPERATURE 

13  TOOTH  SURFACE  FINISH 


0.25000 

30000000.00000 

20.00000 

1.50000 
10.00000 

4500.00000 

0.28300 

9800.00000 

2.50000 
8.93000 
0.90000 

120.00000 

32.00000 


PS  I 

DEGREES 

HORSEPOWER 

RPM 

LBS/IN**3 
PS  I 


DEGREES  F 
RMS 


THERE  ARE  3  INDEPENDENT  DESIGN  VARIABLES. 


ESTIMATED  VALUES: 


LOW  HIGH  INITIAL 


1  PINION  TEETH 

2  DIAMETRAL  PITCH 

3  FACE  WIDTH 


10.00000  100.00000  40.00000 
4.00000  28.00000  10.00000 
0.50000  5.00000  2.50000 


THE  14  CONSTRAINT  LIMITS  ARE: 


1  INVOL.  INTERFERENCE 

2  FACE  WIDTH  TO  DIAMETER 

3  CENTER  DISTANCE 

4  GEAR  AND  PINION  WEIGHT 

5  PINION  TORQUE 

6  TRANSMITTED  FORCE 

7  PITCH  LINE  VELOCITY 

8  TOTAL  DYNAMIC  LOAD 

9  AGMA  BENDING  STRESS 

10  FULL  LOAD  CONTACT  STRESS 

11  GEAR  TIP  HERTZ  PRESSURE 

12  PV  FACTOR 

13  FLASH  TEMPERATURE 

14  MESH  LIFE 


0.00100 

RADIANS 

0.50000 

RATIO 

0.00000 

INCHES 

0.00000 

POUNDS 

0.00000 

LB-IN 

0.00000 

POUNDS 

0.00000 

FT/MIN 

0.00000 

POUNDS 

40000.00000 

PS  I 

150000.00000 

PS  I 

150000.00000 

PS  I 

100.00000 

M  PSI-FT/MIN 

275.00000 

DEGREES  F 

2.00000 

10**3  HOURS 

1.0/INCHES 

INCHES 


TYPE 

LOWER 

UPPER 

LOWER 

LOWER 

LOWER 

LOWER 

LOWER 

LOWER 

UPPER 

UPPER 

UPPER 

UPPER 

UPPER 

LOWER 
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Table  20  Continued 

Spur  Gear  Mesh  Output  File  for  Minimum  Size 


MINIMIZE  THE  OBJECTIVE  FUNCTION,  WHICH  HAS  2  TERMS. 

OBJ  =  THE  LINEAR  SUM  OF: 

1  CENTER  DISTANCE  IN  INCHES  TIMES  1.0000 

2  MESH  LIFE  IN  10**3  HOURS  TIMES  0.0000 


OPTIMIZATION  SUCCESSFUL  IN  33  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 


1  PINION  TEETH 

2  DIAMETRAL  PITCH 

3  FACE  WIDTH 


X(I) 

20.39544 

13.91159  1.0/INCHES 
0.73295  INCHES 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  1.83259 


ITS  COMPONENTS  ARE: 


1  CENTER  DISTANCE  =  1.8326  INCHES  TIMES  1.0000 

2  MESH  LIFE  =  2.0065  10**3  HOURS  TIMES  0.0000 


THE  LAST  CHANGE  IN  THE  OBJECTIVE  FUNCTION  =  -0.326633E-04 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.976563E-04 


THE  14  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1 

INVOL.  INTERFERENCE 

= 

.10861 

RADIANS 

.lOOOOE-02 

LOWER 

2 

FACE  WIDTH  TO  DIAMETER 

= 

.49980 

RATIO 

.50000 

UPPER 

3 

CENTER  DISTANCE 

= 

1.8326 

INCHES 

.00000 

LOWER 

4 

GEAR  AND  PINION  WEIGHT 

1.1376 

POUNDS 

.00000 

LOWER 

5 

PINION  TORQUE 

= 

140.06 

LB-IN 

.00000 

LOWER 

6 

TRANSMITTED  FORCE 

191.07 

POUNDS 

.00000 

LOWER 

7 

PITCH  LINE  VELOCITY 

= 

1727.1 

FT/MIN 

.00000 

LOWER 

8 

TOTAL  DYNAMIC  LOAD 

= 

418.72 

POUNDS 

.00000 

LOWER 

9 

AGMA  BENDING  STRESS 

= 

25097. 

PSI 

40000. 

UPPER 

10 

FULL  LOAD  CONTACT  STRESS 

= 

.14999E+06 

PS  I 

.15000E+06 

UPPER 

11 

GEAR  TIP  HERTZ  PRESSURE 

= 

.12518E+06 

PSI 

.15000E+06 

UPPER 

12 

PV  FACTOR 

= 

92.633 

M  PSI-FT/MIN 

100.00 

UPPER 

13 

FLASH  TEMPERATURE 

= 

216.49 

DEGREES  F 

275.00 

UPPER 

14 

MESH  LIFE 

= 

2.0006 

10**3  HOURS 

2.0000 

LOWER 
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Table  20  Continued 


Sour  Gear  Mesh  Output  File  for  Minimum  Size 


DESIGN  CHECK 


1  PINION  TEETH 

2  DIAMETRAL  PITCH 

3  FACE  WIDTH 


22.00000 

14.00000  1.0/INCHES 
0.75000  INCHES 


THE  MINIMUM  OBJECTIVE  FUNCTION  =  1.96429 


1  CENTER  DISTANCE 

2  MESH  LIFE 


=  1.9643 
=  7.8155 


ITS  COMPONENTS  ARE: 

INCHES  TIMES 

10**3  HOURS  TIMES 


1.0000 

0.0000 


THE  14  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1 

INVOL.  INTERFERENCE 

— 

.12489 

RADIANS 

.lOOOOE-02 

LOWER 

2 

FACE  WIDTH  TO  DIAMETER 

= 

.47727 

RATIO 

.50000 

UPPER 

3 

CENTER  DISTANCE 

= 

1.9643 

INCHES 

.00000 

LOWER 

4 

GEAR  AND  PINION  WEIGHT 

= 

1.3379 

POUNDS 

.00000 

LOWER 

5 

PINION  TORQUE 

= 

140.06 

LB-IN 

.00000 

LOWER 

6 

TRANSMITTED  FORCE 

= 

178.25 

POUNDS 

.00000 

LOWER 

7 

PITCH  LINE  VELOCITY 

= 

1851.3 

FT/MIN 

.00000 

LOWER 

8 

TOTAL  DYNAMIC  LOAD 

= 

398.15 

POUNDS 

.00000 

LOWER 

9 

AGMA  BENDING  STRESS 

= 

22845. 

PS  I 

40000. 

UPPER 

10 

FULL  LOAD  CONTACT  STRESS 

.13873E+06 

PS  I 

.15000E+06 

UPPER 

11 

GEAR  TIP  HERTZ  PRESSURE 

= 

.10980E+06 

PS  I 

.15000E+06 

UPPER 

12 

PV  FACTOR 

= 

81.164 

M  PSI-FT/MIN 

100.00 

UPPER 

13 

FLASH  TEMPERATURE 

= 

206.59 

DEGREES  F 

275.00 

UPPER 

14 

MESH  LIFE 

7.8155 

10**3  HOURS 

2.0000 

LOWER 
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lists  the  found  design  parameter  values,  the  objective  function  values  and  the 
constrained  variable  values  with  their  limits.  Although  it  is  an  "ideal" 
design  with  20.395  teeth  on  the  pinion,  this  theoretical  optimum  identifies 
the  region  of  good  designs.  The  smallest  center  distance  for  a  life  of 
2,000  hours  is  about  one  and  seven-eights  inches  and  the  limiting  constraints 
were:  the  face  width  to  diameter  ratio,  the  full  load  contact  stress  and  the 
mesh  life.  The  watched  variables  in  the  constraint  list  tell  us  the  weight  of 
the  design  and  its  loads  and  velocities. 

The  third  page  is  a  design  check  for  a  design  with  22  teeth  on  the 
pinion  and  a  diametral  pitch  of  14  with  a  face  width  of  0.75  inches.  This 
design  has  a  slightly  larger  center  distance  of  1.964  inches,  but  all 
constraints  are  satisfied  and  the  mesh  life  is  about  7,800  hours.  The  gear 
and  pinion  weigh  1.34  pounds,  the  pinion  torque  is  140  pound-inches,  the 
transmitted  force  is  178  pounds,  and  with  a  pitch  line  velocity  of  1850  feet 
per  minute,  the  dynamic  load  is  estimated  to  be  400  pounds.  Figure  15  is  a 
drawing  of  this  design  showing  the  pinion  and  the  gear  with  its  33  teeth  in 
mesh. 

With  a  few  changes  in  the  input  data  file,  one  can  find  a  design  with 
the  greatest  life  for  a  given  center  distance.  As  can  be  seen  in  the  second 
output  file  of  Table  21,  five  changes  were  made  in  the  input  data  file:  1)  the 
problem  title  was  changed,  2)  the  center  distance  limit  was  changed  from  a 
lower  bound  of  zero  to  an  upper  bound  of  2.5  inches,  3)  the  direction  of 
optimization  was  changed  from  MIN  to  MAX,  and  4  &  5)  the  two  weighting  factors 
in  the  objective  function  were  switched  to  multiply  the  center  distance 
objective  term  by  zero  and  the  mesh  life  objective  term  by  one. 
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MINIMUM  SIZE  SPUR  GEAR  DESIGN 
FIGURE  15 


Table  21 


Sour  Gear  Mesh  Output  File  for  Maximum  Life 

MAXIMUM  LIFE  SPUR  GEAR  WITH  A  FIXED  SIZE 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 


FIXED  DESIGN  REQUIREMENTS: 


1 

POISSON’S  RATIO 

0.25000 

2 

ELASTIC  MODULUS 

30000000.00000 

PSI 

3 

PRESSURE  ANGLE 

20.00000 

DEGREES 

4 

GEAR  RATIO 

1.50000 

5 

TRANSMITTED  POWER 

10.00000 

HORSEPOWER 

6 

PINION  SPEED 

4500.00000 

RPM 

7 

MATERIAL  WEIGHT  DENSITY 

0.28300 

LBS/IN**3 

8 

MATERIAL  SURFACE  CONSTANT 

9800.00000 

PSI 

9 

WEIBULL  SLOPE 

2.50000 

10 

LOAD- LIFE  FACTOR 

8.93000 

11 

RELIABILITY 

0.90000 

12 

BASE  TEMPERATURE 

120.00000 

DEGREES  F 

13 

TOOTH  SURFACE  FINISH 

32.00000 

RMS 

THERE  ARE  3  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 


LOW  HIGH  INITIAL 


1  PINION  TEETH 

10.00000  100.00000  40.00000 

2  DIAMETRAL  PITCH 

4.00000  28.00000  10.00000 

1.0/INCHES 

3  FACE  WIDTH 

0.50000  5.00000  2.50000 

INCHES 

THE 

14  CONSTRAINT  LIMITS  ARE: 

TYPE 

1 

INVOL.  INTERFERENCE 

0.00100 

RADIANS 

LOWER 

2 

FACE  WIDTH  TO  DIAMETER 

0.50000 

RATIO 

UPPER 

3 

CENTER  DISTANCE 

2.50000 

INCHES 

UPPER 

4 

GEAR  AND  PINION  WEIGHT 

0.00000 

POUNDS 

LOWER 

5 

PINION  TORQUE 

0.00000 

LB-IN 

LOWER 

6 

TRANSMITTED  FORCE 

0.00000 

POUNDS 

LOWER 

7 

PITCH  LINE  VELOCITY 

0.00000 

FT/MIN 

LOWER 

8 

TOTAL  DYNAMIC  LOAD 

0.00000 

POUNDS 

LOWER 

9 

AGMA  BENDING  STRESS 

40000.00000 

PSI 

UPPER 

10 

FULL  LOAD  CONTACT  STRESS 

150000.00000 

PSI 

UPPER 

11 

GEAR  TIP  HERTZ  PRESSURE 

150000.00000 

PSI 

UPPER 

12 

PV  FACTOR 

100.00000 

M  PSI-FT/MIN 

UPPER 

13 

FLASH  TEMPERATURE 

275.00000 

DEGREES  F 

UPPER 

14 

MESH  LIFE 

0.00000 

10**3  HOURS 

LOWER 
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Table  21  Continued 


Sour  Gear  Mesh  Output  File  for  Maximum  Life 
MAXIMIZE  THE  OBJECTIVE  FUNCTION,  WHICH  HAS  2  TERMS. 

OBJ  =  THE  LINEAR  SUM  OF: 

1  CENTER  DISTANCE  IN  INCHES  TIMES 

2  MESH  LIFE  IN  10**3  HOURS  TIMES 


OPTIMIZATION  SUCCESSFUL  IN  48  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 

X(I) 

1  PINION  TEETH  26.19858 

2  DIAMETRAL  PITCH  13.09929  1.0/INCHES 

3  FACE  WIDTH  1.00000  INCHES 


THE  MAXIMUM  OBJECTIVE  FUNCTION  =  4689.19  ,  ITS  COMPONENTS  ARE: 

1  CENTER  DISTANCE  =  2.5000  INCHES  TIMES  0.0000 

2  MESH  LIFE  =  4689.2  10**3  HOURS  TIMES  1.0000 


THE  LAST  CHANGE  IN  THE  OBJECTIVE  FUNCTION  =  0.976563E-03 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.953674E-07 


THE  14  CONSTRAINT  VALUES  ARE: 

LIMIT  TYPE 


1 

INVOL.  INTERFERENCE 

= 

.15891 

RADIANS 

.lOOOOE-02 

LOWER 

2 

FACE  WIDTH  TO  DIAMETER 

= 

.50000 

RATIO 

.50000 

UPPER 

3 

CENTER  DISTANCE 

2.5000 

INCHES 

2.5000 

UPPER 

4 

GEAR  AND  PINION  WEIGHT 

= 

2.8895 

POUNDS 

.00000 

LOWER 

5 

PINION  TORQUE 

= 

140.06 

LB- IN 

.00000 

LOWER 

6 

TRANSMITTED  FORCE 

=: 

140.06 

POUNDS 

.00000 

LOWER 

7 

PITCH  LINE  VELOCITY 

z= 

2356.2 

FT/MIN 

.00000 

LOWER 

8 

TOTAL  DYNAMIC  LOAD 

= 

334.97 

POUNDS 

.00000 

LOWER 

9 

AGMA  BENDING  STRESS 

12766. 

PS  I 

40000. 

UPPER 

10 

FULL  LOAD  CONTACT  STRESS 

= 

96590. 

PS  I 

.15000E+06 

UPPER 

11 

GEAR  TIP  HERTZ  PRESSURE 

= 

70072. 

PS  I 

.15000E+06 

UPPER 

12 

PV  FACTOR 

=: 

56.012 

M  PSI-FT/MIN 

100.00 

UPPER 

13 

FLASH  TEMPERATURE 

175.53 

DEGREES  F 

275.00 

UPPER 

14 

MESH  LIFE 

= 

4689.2 

10**3  HOURS 

.00000 

LOWER 

0.0000 
I. 0000 
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Table  21  Continued 


DESIGN  CHECK 

X(I) 

1  PINION  TEETH  24.00000 

2  DIAMETRAL  PITCH  12.00000  1.0/INCHES 

3  FACE  WIDTH  1.00000  INCHES 


THE  MAXIMUM  OBJECTIVE  FUNCTION  =  4432.67  ,  ITS  COMPONENTS  ARE: 

1  CENTER  DISTANCE  =  2.5000  INCHES  TIMES  0.0000 

2  MESH  LIFE  =  4432.7  10**3  HOURS  TIMES  1.0000 


THE  14  CONSTRAINT  VALUES  ARE: 


LIMIT 

TYPE 

1 

INVOL.  INTERFERENCE 

.14243 

RADIANS 

.lOOOOE-02 

LOWER 

2 

FACE  WIDTH  TO  DIAMETER 

=: 

.50000 

RATIO 

.50000 

UPPER 

3 

CENTER  DISTANCE 

= 

2.5000 

INCHES 

2.5000 

UPPER 

4 

GEAR  AND  PINION  WEIGHT 

=r 

2.8895 

POUNDS 

.00000 

LOWER 

5 

PINION  TORQUE 

140.06 

LB-IN 

.00000 

LOWER 

6 

TRANSMITTED  FORCE 

140.06 

POUNDS 

.00000 

LOWER 

7 

PITCH  LINE  VELOCITY 

= 

2356.2 

FT/M IN 

.00000 

LOWER 

8 

TOTAL  DYNAMIC  LOAD 

= 

334.97 

POUNDS 

.00000 

LOWER 

9 

AGMA  BENDING  STRESS 

r: 

12009. 

PS  I 

40000. 

UPPER 

10 

FULL  LOAD  CONTACT  STRESS 

= 

97085. 

PS  I 

.15000E+06 

UPPER 

11 

GEAR  TIP  HERTZ  PRESSURE 

== 

73216. 

PS  I 

.15000E+06 

UPPER 

12 

PV  FACTOR 

= 

63.518 

M  PSI-FT/MIN 

100.00 

UPPER 

13 

FLASH  TEMPERATURE 

= 

180.12 

DEGREES  F 

275.00 

UPPER 

14 

MESH  LIFE 

= 

4432.7 

10**3  HOURS 

.00000 

LOWER 
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For  this  new  objective  of  maximizing  the  life  of  a  larger  gear  set,  the 
program  found  an  optimum  design  with  a  life  of  4.69  million  hours  in  48  steps. 
This  design  was  bounded  by  the  center  distance  limit  of  2.5  inches  and  the 
length  to  diameter  ratio  limit  of  0.5.  This  optimum  also  had  an  unrealistic 
number  of  teeth  on  the  pinion  of  26.2  and  diametral  pitch  of  13.1.  The  last 
page  of  Table  21  shows  the  results  of  a  nearby  design  check  using  24  teeth  on 
the  pinion,  a  diametral  pitch  of  12  and  a  face  width  of  one  inch.  This  design 
is  shown  in  Figure  16  with  24  teeth  on  the  pinion  and  36  teeth  on  the  gear  and 
the  requested  2.5  inch  center  distance. 

The  realistic  design  has  a  slightly  lower  life  of  4.43  million  hours,  a 
gear  and  pinion  weight  of  2.9  pounds,  a  pinion  torque  of  140  pound-inches,  a 
transmitted  force  of  140  pounds,  a  pitch  line  velocity  of  2356  feet  per 
minute,  and  an  estimated  dynamic  load  of  335  pounds.  Since  it  is  larger  than 
the  minimum  size  design,  it  also  has  lower  loads  and  stresses. 
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MAXIMUM  LIFE  SPUR  GEAR  DESIGN 
FIGURE  16 


106 


WEIBULL  DATA  FITTING 


The  fourth  example  application  is  the  use  of  the  optimizer  to  fit 
experimental  life  data  to  a  two-parameter  Weibull  distribution.  With  its 
abilities  to  iterate  in  n  independent  variables  and  provide  an  easy  check  of 
the  performance  of  alternate  solutions,  the  computer  program  SEEK. FOR  offers  a 
variety  of  uses.  One  valuable  application  is  that  of  fitting  a  curve  to 
experimental  data. 

Theory 

A  nonlinear  least  squares  fitting  procedure  compares  the  measured  data 
to  the  fitted  curve  at  each  point.  The  sum  of  the  squares  of  the  errors 
between  the  measured  data  and  the  values  on  the  fitted  curve  is  a  positive 
definite  scalar  measure  of  the  scatter  of  the  data  about  the  curve.  Taking 
the  square  root  of  this  squared  error  sum  divided  by  the  number  of  data  points 
less  one  gives  a  dispersion  with  the  same  dimension  as  the  measured  quantity 
[22].  For  the  two-parameter  Weibull  relationship  of  equation  (33),  this 
dispersion  would  be: 


(^c  - 

Nq-  1 


1/2 


(38) 


where  is  the  reliability  on  the  curve,  is  the  median  rank  reliability  of 
the  measured  life  and  Nq  is  the  number  of  data  points. 

Thus  one  can  use  SEEK. FOR  to  determine  the  ninety-percent  reliability 
life,  fjQ,  and  Weibull  slope,  b,  from  a  set  of  life  test  data.  Table  22 
summarizes  this  problem  in  the  design  optimization  format  with:  the  lone 
constant  being  the  number  of  data  points;  the  two  design  parameters,  and 
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Table  22 

Sour  Gear  Mesh  Optimization  Parameters 


Constants 

Design  Parameters 

Inequality  Constraints 

Objective  Function 

*^0 

^10 

b 

Np  >  0 

^^D^min 

Table  23 


Median  Ranked  Life  Data 
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b;  the  inequality  constraint  of  having  a  positive  number  of  data  points;  and 
the  objective  function  of  minimizing  the  dispersion, 

For  a  group  of  identical  units  subjected  to  an  identical  life  test, 
different  units  will  fail  at  different  times.  By  recording  the  times  to 
failure  and  ranking  them  from  the  shortest  to  the  longest,  one  can  use  a 
median  rank  table  [23]  to  list  the  median  reliability  and  test  life  in  a  table 
or  data  file  similar  to  the  numerical  values  in  Table  23.  This  table  lists 
the  time  to  failure  and  the  corresponding  median  reliability  for  a  life  test 
of  fifteen  units. 

Programming 

Inside  SEEK. FOR  is  the  additional  capability  to  read  a  second  data  file 
in  addition  to  the  input  data  file.  This  action  is  enabled  by  setting  the 
first  constant  in  the  constant  array  to  the  number,  n,  of  data  point  pairs, 
(XD^,YD^),  in  the  data  file  and  including  the  words,  "DATA  POINTS,"  in  its 
description.  When  this  is  combined  with  the  presence  of  a  data  file  of  n 
lines  with  n  point  pairs;  XD^  ,  YD^;  and  a  name  which  has  the  same  prefix  as 
the  ".IN"  file  but  the  extension  ".DAT;"  SEEK  will  open  and  read  the  ".DAT" 
file  after  it  has  read  the  ".IN"  file.  The  data  pairs  will  be  placed  in  a 
common  block,  COMMON/CURVE/XDP(200) ,YDP(200) ,  for  use  in  subroutine  VALUES. 

If  the  limit  of  200  data  pairs  is  too  low,  changing  this  dimension  in  the 
program  will  allow  larger  data  sets  to  be  fit. 

Table  24  lists  subroutines  BOUNDS  and  VALUES  which  will  fit  a  two- 
parameter  Weibull  distribution  to  the  data  in  the  ".DAT"  file.  Subroutine 
BOUNDS  checks  that  the  number  of  data  points  is  positive,  since  it  has  to  do 
something  to  permit  the  overall  program  to  execute.  Subroutine  VALUES  uses 
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Table  24 


Two-Parameter  Weibull  Fitting  Analysis  Subroutines 
SUBROUT I NE  BOUNDS ( CONST , NCO , X , NX , VCSTR , NCS ) 

INPUT 

CONST(l)  =  THE  NUMBER  OF  DATA  POINTS 
OUTPUT 

VCSTR(l)  =  THE  NUMBER  OF  DATA  POINTS 

DIMENSION  X(NX) ,CONST(NCO) ,VCSTR(NCS) 

VCSTR(I)  =  CONST(l) 

RETURN 

END 

SUBROUTINE  VALUES ( CONST , NCO , X , NX , OBJECT , NOB) 

MINIMIZE  THE  SUM  OF  THE  FITTED  ERRORS  SQUARED 
INPUTS 

CONST(I)  =  THE  NUMBER  OF  DATA  POINTS 

X(l)  =  RL  -  THE  SOUGHT  90%  RELIABILITY  LIFE 

OF  THE  DISTRIBUTION 
X(2)  =  B  -  THE  SOUGHT  WEIBULL  SLOPE 

OUTPUT 

OBJECT(l)  =  THE  DEVIATION  IN  THE  ERRORS 

DIMENSION  X(NX) ,CONST(NCO) ,OBJECT(NOB) 
COMMON/CURVE/TDP(200) , RDP(200) 

N  =  CONST(l) 

RL  =  X(l) 

B  =  X(2) 

ERS  =  0.0 
DO  10  I  =  I,N 

R  =  I.0/EXP(AL0G(I.0/0.9)*(TDP(I)/RL)**B) 

E  =  R  -  RDP(I) 

ERS  =  ERS  +  E*E 
10  CONTINUE 
XN  =  N  -  1 

OBJECT(l)  =  SQRT(ERS/XN) 

RETURN 

END 
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equations  (33)  and  (38)  to  determine  the  dispersion  of  the  measured  data  about 
the  trial  distribution. 

Numerical  Example 

To  illustrate  this  data  fitting  process,  the  fifteen  data  points  of 
Table  23  and  the  analysis  subroutines  of  Table  24  will  be  used.  A  listing  of 
the  input  data  file  "W.IN"  which  will  initiate  the  curve  fitting  is  given  in 
Table  25.  In  this  file;  the  first  line  is  the  problem  title;  the  second  line 
is  the  number  of  constants,  which  is  one;  and  the  next  three  lines  are  the 
number  of  data  points,  15;  the  description  "NUMBER  OF  DATA  POINTS"  which 
includes  the  words  "DATA  POINTS"  and  a  blank  line  for  the  units.  Following 
this  is  a  line  with  the  number  of  independent  parameters,  two.  The  next  six 
lines  give  the  two  distribution  parameters,  with  their  low,  high  and  initial 
estimates.  This  is  followed  by  the  number  of  design  constraints,  one,  and 
three  lines  which  describe  the  lower  bound  on  the  number  of  data  points.  On 
the  next  line  are  the  three  letters  "MIN"  to  signal  that  the  objective 
function  is  to  be  minimized  and  four  lines  which  give  the  number  of  object 
functions  as  one  and  the  weighting  coefficient,  description  and  units  for  the 
reliability  error  dispersion.  Once  again,  the  unit  line  is  blank. 

Compiling  the  subroutines  of  Table  24,  linking  them  with  SEEK. OBJ,  and 
running  the  resulting  program  with  the  two  data  files  W.DAT  which  holds  the 
numbers  of  Table  23  and  W.IN  which  is  listed  in  Table  25  yields  the  output 
file  of  Table  26.  This  file  echoes  the  input  file  information  and  reports  a 
successful  curve  fit  in  43  steps.  The  fitted  curve  has  an  life  of 
238.8  thousand  hours  and  a  Weibull  slope  of  2.296  with  a  dispersion 
reliability  of  0.0317. 


Ill 


Table  25 


Two-Parameter  Wei  bull  Fitting  Input  Data  File 


WEIBULL  RELIABILITY  DISTRIBUTION  FIT 
1 

15 

NUMBER  OF  DATA  POINTS 
2 

100.0  2000.0  500.0 
90%  RELIABILITY  LIFE 
10**3  HOURS 
0.0  6.0  2.0 
WEIBULL  SLOPE 

1 

LOWER  0.0 

NO.  OF  DATA  POINTS 

MIN 

1 

1.0 

RELIABILITY  ERROR  DISPERSION 
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WEIBULL  RELIABILITY  DISTRIBUTION  FIT 

DESIGN  WITH  MODIFIED  GRADIENT  OPTIMIZATION 
USING  A  MAXIMUM  STEP  LIMIT  AND  SCALED  VARIABLES. 

FIXED  DESIGN  REQUIREMENTS: 

1  NUMBER  OF  DATA  POINTS  15.00000 

THERE  ARE  2  INDEPENDENT  DESIGN  VARIABLES. 

ESTIMATED  VALUES: 

LOW  HIGH  INITIAL 

1  90%  RELIABILITY  LIFE  100.0000  2000.0000  500.0000  10**3  HOURS 

2  WEIBULL  SLOPE  0.0000  6.0000  2.0000 

THE  1  CONSTRAINT  LIMITS  ARE: 

TYPE 

1  NO.  OF  DATA  POINTS  0.00000  LOWER 

MINIMIZE  THE  OBJECTIVE  FUNCTION. 

OBJ  =  RELIABILITY  ERROR  DISPERSION  IN  TIMES  1.0000 

15  POINTS  READ  FROM  TABULAR  DATA  FILE 

OPTIMIZATION  SUCCESSFUL  IN  43  STEPS 
THE  FINAL  DESIGN  VECTOR  IS: 

X(I) 

1  90%  RELIABILITY  LIFE  238.80429  10**3  HOURS 

2  WEIBULL  SLOPE  2.29637 

THE  MINIMUM  OBJECTIVE  FUNCTION  =  0.316652E-01  ,  ITS  COMPONENTS  ARE: 

1  RELIABILITY  ERROR  DISPERSION  =0.31665E-01  TIMES  1.0000 

THE  LAST  CHANGE  IN  THE  OBJECTIVE  FUNCTION  =  -0.372529E-08 

THE  LAST  STEP  CHANGE  SIZE  FOR  THE  DESIGN  VARIABLE  =  0.610352E-05 


Table  26  Continued 

Two-Parameter  Wei bu11  Fitting  Output  Data  File 


THE  1  CONSTRAINT  VALUES  ARE: 
1  NO.  OF  DATA  POINTS  =  15.000 


LIMIT  TYPE 

.00000  LOWER 


DESIGN  CHECK 


X(I) 


90%  RELIABILITY  LIFE 

WEIBULL  SLOPE 

240.00000 

2.30000 

10**3  HOURS 

MINIMUM  OBJECTIVE  FUNCTION  = 

0.315955E-01  , 

ITS  COMPONENTS  ARE 

• 

RELIABILITY  ERROR  DISPERSION 

=0.31596E-01 

TIMES 

1.0000 

THE  1  CONSTRAINT  VALUES  ARE: 

LIMIT 

TYPE 

1.  OF  DATA  POINTS  =  15. 

000 

.00000 

LOWER 

II4 


As  shown  in  the  design  check,  rounding  the  fjQ  life  to  240  thousand 
hours  and  the  WeibuH  slope  to  2.3  changes  this  dispersion  to  0.0316.  This 
rounded  function  fits  the  data  slightly  better  than  the  "ideal"  function  as 
found  by  the  optimizer  due  to  the  finite  step  size  procedure  of  the  gradient 
algorithm.  Selecting  a  different  starting  position  or  sensitivity  ranges  for 
the  two  Weibull  parameters  could  result  in  an  "ideal"  solution  which  is 
slightly  better  than  the  rounded  solution,  but  it  would  not  be  significantly 
different.  A  good  fit  is  achieved  with  two  digit  precision  parameters,  using 
the  combination  of  the  optimizer  and  the  design  check. 

Figure  17  is  a  two-parameter  Weibull  graph  of  the  measured  data  with  its 
median  rank  reliability  plotted  as  crosses.  On  this  plot  is  drawn  the  rounded 
two-parameter  Weibull  distribution  which  has  been  fit  to  the  data.  It  is 
always  a  good  idea  to  draw  the  fitted  curve  with  the  data  points  to  which  it 
is  fit.  Since  the  dispersion  is  a  single  value  representing  the  overall  fit 
of  the  curve,  some  local  anomalies  may  exist  in  a  fitted  curve.  A  graph  shows 
these.  If  the  data  does  not  have  the  general  shape  of  the  fitted  function, 
either  a  partial  fit  in  the  region  of  interest  or  a  more  sophisticated 
function  should  be  used. 

By  changing  the  ranking  from  median  to  both  low  and  high  reliability 
ranks,  the  program  can  be  run  twice  more  to  determine  confidence  bounds  on  the 
distribution.  These  two  runs  would  give  the  user  a  low  reliability  fit  to  the 
test  data  and  a  high  reliability  fit  to  the  same  data.  Plotting  these  two 
curves  on  the  graph  of  Figure  17  would  show  the  confidence  range  for  the 
distribution. 

To  investigate  the  potential  improvement  in  fitting  the  data  with  a 
higher  order  Weibull  distribution,  a  three-parameter  Weibull  distribution 
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.0  200.0  500.0  1,000.0 

LIFE-  I  (10^  HOURS) 

TWO-PARAMETER  WEIBULL  DISTRIBUTION 
WITH  Q=  240  *  1 0^  HOURS  AND  b  =  2.3 


FIGURE  17 


could  be  fit  to  the  same  data  by  changing  the  reliability  equation  and  adding 


a  third  independent  parameter  to  the  optimization.  The  equation  for  the  three 
parameter  Wei  bull  distribution  is; 


Ln 


1 

R 


=  Ln 


1 

0.9 


1. 


t. 


(39) 


*10  '0 

where  the  minimum  life,  is  the  additional  independent  parameter.  The 
third  parameter,  replaces  zero  as  the  minimum  life  that  all  units  would 
realize. 

Modifying  the  program  and  the  input  data  file  would  enable  one  to  obtain 
a  three-parameter  Wei  bull  distribution  fit  to  the  same  data.  Comparing  the 
two  dispersions  would  indicate  whether  the  additional  complexity  of  the  three- 


parameter  distribution  is  justified  in  the  life  model.  A  small  minimum  life 
would  confirm  the  adequacy  of  the  two-parameter  Weibull  distribution  as  well. 
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OPTIMIZATION  METHOD 


Parameter  Scaling 

As  illustrated  in  the  examples,  the  optimization  procedure  begins  with 
several  vectors.  An  optimization  solution  is  the  design  variable  vector,  X, 
which  minimizes  or  maximizes  the  objective  function  value,  M,  with  all 
constraint  values,  Vj^,  bounded  by  their  specified  limits.  A  procedure  starts 
with  a  guess  for  the  design  variables,  X,  and  iterates  using  gradients  to  find 
the  optimal  values.  Opportunity  is  then  given  to  the  user  to  try  alternate 
solutions  and  compare  their  properties. 

To  maintain  balance  among  the  independent  design  parameters,  the  design 
space  is  scaled  into  a  dimensionless  design  space  [3].  The  scaled  design 
parameters,  Y^,  vary  from  -  1.0  to  +  1.0  as  specified  by  upper  and  lower 
bounds  on  the  independent  design  parameters,  X^,  such  that: 


-  1.0  <  Y.  <  +  1.0 


as 


X 


Li 


< 


*1 


< 


X 


Ui 


(40) 

(41) 


This  linear  transformation,  which  is  shown  in  Figure  18,  is: 


where 


Yi  =  d.  X.  .  b. 


d,  = 


^Ui  “  ^Li 


(42) 


(43) 


and 


Xui  ^  X,. 
^Ui  ■  ^Li 


(44) 
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The  actual  design  variable,  X.,  can  be  retrieved  from  the  scaled  variable,  Y., 
by  the  inverse  transformation: 

Y.  -  b. 

X  =  !j _ I  (45) 


Gradients 

Central  to  this  method  is  the  gradient  calculation  which  is  performed 
with  small  perturbations  in  the  design  variables  from  the  nominal  position. 
The  gradient  in  the  merit  function,  V  M,  is  calculated  [1]  as: 


V  M 


aM 

aM 

ar^ 

aM 


(46) 


where. 


aM 

3Y7 


M(Yj,-,Y.+AY,-,Y^)  -  M(Y] 


,Y.,-,Y^) 


'W 


(47) 


In  the  program,  the  small  change  AY,  which  is  made  in  each  Y^,  is  set  at  0.001 
which  is  0.05  percent  of  the  full  range  of  a  scaled  design  parameter. 

The  magnitude  of  the  gradient  vector  is  given  by: 
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n 


2  1/2 


V  M 


(48) 


In  the  program,  the  step  size  is  separated  from  the  gradient  magnitude  but  for 
minimization  its  direction  is  opposite  to  the  gradient  direction,  so  negative 
unit  gradient  vectors  are  determined: 


I  VM  I 

For  maximization,  the  sign  in  equation  (49)  reverses.  This  sign  reversal 
also  occurs  for  the  constraint  gradient  when  the  limit  is  a  lower  bound. 
Search  Directions 

In  the  simple  gradient  method  which  is  used  in  the  acceptable  design 
region,  equation  (49)  defines  the  direction  of  change  in  the  scaled  design 
vector,  Y.. 


Y .  +  AS  ■  V  m 


(50) 


If  no  constraints  are  violated,  this  will  be  the  next  value  for  Y^  in  the 
search. 

For  stability  and  directness,  a  nominally  fixed  step  size,  AS,  is  used 
in  this  optimization.  Initially,  the  step  size  is  0.1,  which  is  five  percent 
of  the  range  of  a  single  design  parameter.  Then,  whenever  a  local  minimum  is 
reached  or  the  search  is  trapped  in  a  constraint  corner,  the  procedure  halves 
the  step.  To  complete  the  search  procedure,  the  program  declares  a  solution 
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when  the  percent  change  in  the  merit  function  is  less  than  a  pre  set  limit  of 


0.0001: 


<  0.0001 


(51) 


Two  other  modes  of  searching  are  used  in  the  procedure:  1)  poor  guess 


correction  when  constraints  are  violated  initially,  and  2)  feasible  direction 
searching  near  constraint  boundaries. 

These  modes  are  enabled  with  a  second  gradient.  Just  as  one  can 
calculate  the  gradient  in  the  objective  function,  one  also  can  calculate  the 
gradient  in  a  violated  constraint  variable: 


I  VVk 


(52) 


For  upper  bound  constraints,  moving  through  the  design  space  in  the  direction 
of  V  V.  will  reduce  the  constraint  value  For  lower  bound  constraints,  a 

IN 

sign  reversal  in  equation  (52)  produces  an  increase  in  the  constraint  value, 
y.,  for  motion  in  the  gradient  direction.  The  vector  sum  of  the  gradients  in 

K. 

the  violated  constraints,  V  h,  is  the  second  gradient  of  the  algorithm: 


E  ^''k 

k 


The  gradient  in  the  violated  constraints,  V  h,  points  towards  the 
acceptable  design  space  from  the  unacceptable  design  space.  By  itself,  it 
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enables  the  algorithm  to  turn  an  unacceptable  initial  guess  into  an  acceptable 
design  trial  by  a  succession  of  steps: 


+  AS  •  V  h 


(54) 


Once  inside  the  acceptable  design  region,  the  algorithm  proceeds  along 
the  steepest  descent  direction  with  equation  (50)  until  the  calculated  step 
places  the  next  trial  outside  the  acceptable  design  space.  To  avoid  this 
condition,  the  algorithm  selects  a  feasible  direction  for  the  next  step. 

Figure  19  shows  a  sloped  constraint  intersecting  vertical  contour  lines  of  the 
objective  function.  This  figure  is  an  enlargement  of  a  small  region  in 
Figure  2,  which  is  a  plot  of  the  length  versus  diameter  for  a  bushing.  In 
this  example,  the  objective  function  is  directly  proportional  to  the  bushing 
diameter  and  independent  of  its  length.  So  a  steepest  descent  direction  is 
always  horizontal  for  the  problem.  The  sloped  constraint  curve  is  the  length 
to  diameter  upper  bound.  Figure  19  shows  vertical  contour  lines  in  the 
frictional  torque  objective  function  and  unit  gradient  vectors  in  this 
objective  function,  V  m,  and  the  impending  constraint,  V  h.  The  two  gradient 
vectors  are  added  at  the  last  viable  design  step.  The  feasible  direction 
selected,  V  f,  is  the  unit  vector  sum  of  these  two  gradients: 


V  m  +  V  h 

V  f  =  - 

V  m  +  V  h 


(55) 


And  the  next  step  becomes: 

Vl  =  Y.  .  AS-Vf 


(56) 
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FEASIBLE  DIRECTION  GRADIENT 


FIGURE  19 


Program  Structure 

In  writing  the  optimization  code,  subroutines  were  used  to  modularize 
the  programming  and  perform  the  vector  calculations  in  a  structured  way. 

Table  1,  which  is  cited  in  the  PROGRAMMING  section  of  this  report,  lists  the 
subroutines  which  compose  this  program  with  short  descriptions  of  their 
functions.  Each  subroutine’s  name  is  related  to  its  operation  to  describe  its 
use. 

The  scaling  and  unsealing  of  the  design  parameter  vector  are  performed 
by  two  routines:  SIZE  and  RESIZE,  so  that  this  linear  transformation  can  be 
performed  anywhere  in  the  program  with  relative  ease.  Subroutines  GRADNT  and 
UNIT  evaluate  and  normalize  any  gradient  vector.  Subroutine  BOUNCE  finds  the 
gradient  sum  of  the  violated  constraints  at  any  design  position.  Subroutine 
CHECK  compares  a  constraint  value  to  its  limit  and  subroutine  WALL  evaluates  a 
specific  constraint  in  the  same  way  that  subroutine  MERIT  evaluates  the 
objective  function. 

At  a  higher  level  in  the  program  structure,  subroutine  BOUNCE  directs 
the  search  for  the  acceptable  design  space  when  the  initial  design  guess 
violates  at  least  one  constraint.  Subroutine  SCAN  performs  the  search  for 
better  designs  within  the  acceptable  design  space  and  subroutine  SCOUT  checks 
the  next  potential  design  for  constraint  violations.  Finally,  subroutines 
BOUNDS  and  VALUES  are  the  user  supplied,  problem  specific  routines  which 
evaluate  the  design  constraints  and  the  objective  function  values  for  the 
design. 

Figure  20  is  a  logic  flow  chart  for  the  optimization  program  which 
includes:  1)  the  reading  of  the  input  data  files,  2)  the  echoing  of  the  input 
data  file  information  at  the  start  of  the  output  data  file,  3)  the  initial 
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READ  DATA  Sc  SET  CONSTANTS 
WRITE  INPUT  DATA  TO  OUTPUT  FILE 


FOR  FIRST  CONSTANT  OF  "DATA  POINTS" 
READ  TABLE  OF  XDP,  YDP  VALUES 


CHECK  CONSTRAINT  Sc  MERIT  FUNCTION  VALUES 


ARE  ALL  CONSTRAINTS  SATISFIED  ? 


CALCULATE  Vm  Sc  LOCATION  OF  NEXT  POINT 
ALONG  Vm  Sc  CONSTRAINT  VALUE  THERE 


ARE  ALL  CONSTRAINTS  SATISFIED  ? 


DID  MERIT  FUNCTION  BECOME  ZERO  ? 
IS  ITS  PERCENT  CHANGE  SMALL  ? 
OR  HAVE  500  TRIALS  BEEN  MADE  ? 


SEEK  PROGRAM  FLOW  CHART 


FIGURE  20 
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search  for  a  good  starting  design,  4)  the  main  optimization  loop,  and  5)  the 
design  verification  loop. 

The  search  for  a  good  starting  design  which  satisfies  all  constraints  is 
performed  in  a  DO  loop  with  a  limit  of  twenty  corrections.  If  a  valid  design 
is  not  found  in  twenty  trials,  both  the  first  and  twentieth  trials  are  written 
to  the  output  file  with  an  analysis  of  all  constraints,  the  file  is  closed  and 
the  program  is  stopped. 

Once  a  good  starting  design  is  found,  the  search  for  the  optimal  design 
is  conducted  in  the  main  DO  loop  with  a  limit  of  five-hundred  iterations.  In 
this  DO  loop,  subroutine  SCOUT  checks  the  constraint  values  at  the  next  design 
point  using  the  objective  function  gradient  increment  as  shown  in 
equation  (50)  to  locate  the  next  design  point.  If  all  the  constraints  are 
satisfied,  the  full  design  step  is  taken  in  the  direction  of  the  objective 
function  gradient  using  equation  (50).  If  at  least  one  constraint  is  not 
satisfied,  subroutine  SCAN  increments  the  design  along  the  feasible  direction 
gradient  using  equation  (56).  If  all  the  constraints  are  still  not  satisfied, 
the  step  size  is  reduced  and  the  half  step  is  taken  in  the  feasible  direction. 

After  a  step  is  taken,  the  objective  function  value  is  then  checked.  If 
it  increased,  the  step  size  is  divided  by  two.  Otherwise,  no  change  is  made. 
Finally,  the  percent  change  in  the  objective  function  is  checked  using  MERIT. 
If  this  change  is  below  the  desired  limit,  a  solution  is  declared  and  the 
procedure  leaves  the  DO  loop.  If  the  change  is  greater  than  the  desired 
limit,  the  DO  loop  indexes  and  the  iteration  process  is  repeated. 

If  the  merit  function  criteria  are  not  satisfied  after  five-hundred 
iterations,  the  program  writes  the  five-hundredth  design  to  the  screen  and  the 
output  file  with  an  error  message  and  provides  the  user  with  the  opportunity 
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to  change  the  design  and  have  it  re-analyzed.  This  documentation  and 
opportunity  for  modification  are  also  given  to  the  user  at  the  end  of  a 
successful  design  optimization  search.  The  final  design  modification  loop  is 
controlled  by  the  user  with  no  termination  count.  Additional  designs  may  be 
evaluated  until  the  user  chooses  to  end  the  program. 
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DISCUSSION  OF  RESULTS 


In  this  report,  a  computer  program,  ’SEEK. FOR,’  is  described  and  its  use 
is  illustrated  with  four  different  examples.  SEEK  has  been  written  to  make 
optimization  available  to  the  general  technical  community  with  only  the  need 
for  modeling  the  applied  problem.  The  program  is  a  general  purpose  optimizing 
tool  which  can  find  the  best  combination  of  design  parameter  values  which 
satisfy  a  series  of  constraints  placed  on  calculable  properties  of  the  design 
parameters.  The  set  of  design  parameter  values  is  ’best’  in  that  it  maximizes 
or  minimizes  a  linear  sum  of  objective  function  terms. 

Although  the  computer  program  requires  the  user  to  write  two  analysis 
subroutines  and  an  input  data  file  for  his  or  her  application,  it  is  written 
to  be  interactive  and  to  communicate  clearly  to  the  user.  By  requiring  text 
labels  for  all  variables,  the  program  is  able  to  label  its  restatement  of  the 
problem  and  all  results  in  the  words  of  the  user.  In  the  input  file,  the  user 
is  able  to  change  the  starting  design  values  and  the  relative  sensitivity  of 
the  individual  design  parameters  by  changing  the  low  and  high  estimates  of 
these  values. 

Based  on  a  bounded  step  gradient  method,  the  program  searches  with 
gradients  in  three  functions:  1)  the  objective  function,  2)  the  violated  or 
nearly  violated  constraints  and  3)  a  vector  sum  of  the  first  two.  The  first 
gradient  provides  a  rapid  solution  path  when  no  obstacles  are  present.  The 
second  gradient  permits  the  user  to  start  with  an  initial  trial  design  which 
may  not  satisfy  all  constraints.  And  the  third  gradient  allows  the  algorithm 
to  continue  improving  the  objective  function  even  when  the  direct  path  is 
blocked  with  one  or  more  constraints. 
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To  enable  the  gradient  calculations,  the  provided  models  for  the 
objective  function  terms  and  the  constrained  variables  need  to  be  continuous. 
Discrete  values  can  be  entered  later  in  a  design  check  mode  which  the  program 
enters  once  a  numerical  optimum  has  been  found  in  the  continuous  design  space. 
This  design  check  mode  is  part  of  the  interactive  operation  of  the  program. 

By  allowing  the  user  to  try  alternate  designs  after  the  numerical  optimum  has 
been  found,  the  program  assists  the  user  in  finding  practical  optimums  which 
satisfy  discrete  requirements  which  can  not  be  treated  in  the  continuous 
models.  An  output  file  is  written  during  the  session  to  document  the  design 
problem  and  all  solutions  which  are  displayed  on  the  screen.  In  this  file, 
each  solution  has  its  design  parameter  values,  objective  function  values,  and 
constrained  variable  values  and  limits  listed  with  the  users  descriptions  as 
labels. 

This  report  illustrates  the  model  preparation  and  programming  required 
to  use  ’SEEK. FOR’  for  four  examples:  1)  a  bushing  design,  2)  a  spring  design, 
3)  a  gear  design,  and  4)  a  curve  fit.  The  bushing  design  problem  illustrates 
the  determination  of  two  parameters:  the  bushing  length  and  its  diameter 
subject  to  three  inequality  constraints,  a  single  objective  function  and 
discrete  size  requirements.  The  example  also  shows  the  program’s  ability  to 
overcome  a  poor  initial  trial  design. 

The  spring  design  problem  illustrates  the  conversion  of  equality 
constraints  to  simplify  the  problem  by  taking  a  four-parameter  problem  and 
treating  it  as  a  problem  with  two  independent  parameters.  The  number  of 
active  coils  in  the  spring  and  the  free  length  of  the  spring  are  shown  to  be 
dependent  parameters  and  the  wire  diameter  and  mean  coil  diameter  are  left  as 
independent  parameters.  The  example  finds  three  separate  designs  which  all 
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have  the  same  stiffness  and  design  factor  under  the  same  loading  but  which 
satisfy  three  separate  objective  functions.  By  changing  the  objective 
function  for  the  same  design  constraints,  the  program  finds  the  lightest 
spring,  the  smallest  volume  spring  and  the  shortest  spring  in  separate  runs. 

A  visual  comparison  of  the  three  springs  shows  their  obvious  differences. 

The  gear  design  problem  illustrates  a  three  parameter  design  problem 
which  reguires  considerably  more  analysis  than  the  first  two.  The  three 
parameters  are  the  number  of  teeth  on  the  pinion,  the  diametral  pitch  and  the 
face  width  of  the  gears.  Designs  are  obtained  for  a  minimum  size  at  a  given 
life  and  for  a  maximum  life  at  a  given  size.  The  example  illustrates  the  use 
of  additional  subroutines  by  the  two  analysis  routines:  BOUNDS  and  VALUES. 

The  curve  fit  problem  illustrates  the  use  of  the  optimizing  program  as  a 
tool  for  minimizing  the  dispersion  of  the  data  from  a  modeled  function. 
Although  the  program  can  minimize  the  overall  error  between  a  fitted  curve  and 
the  experimental  data,  it  is  recommended  that  a  plot  of  the  fitted  curve  and 
the  original  data  be  made  to  verify  that  the  function  truly  models  the  data 
closely  in  the  region  of  highest  interest.  One  feature  of  this  use  of  the 
program  is  the  ability  in  the  design  check  mode  to  quickly  see  the  influence 
of  curve  coefficients  on  the  goodness  of  fit. 

No  optimization  program  can  solve  all  problems.  The  use  of  this  program 
is  limited  to  solving  problems  which  can  be  modeled  continuously  with  finite 
constraint  and  objective  function  values  over  the  search  area.  The  program 
contains  arbitrary  size  limits  of  forty  for  the  constants  and  constraints  and 
fifteen  for  the  design  parameters  and  objective  function  terms.  As  the  number 
of  independent  design  parameters  increases  from  the  low  numbers  of  these 
examples,  gradient  searching  becomes  more  difficult.  More  steps  are  required 
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to  find  a  solution  and  more  calculations  are  required  for  each  step.  Design 
parameters,  which  do  not  change,  can  dilute  the  effectiveness  of  the  method 
with  additional  calculations.  Often,  a  large  optimization  problem  can  be 
reduced  to  a  simpler  one  which  permits  a  more  rapid  and  direct  solution  by  the 
removal  of  design  parameters  which  do  not  change  from  the  design  parameter 
vector. 

For  any  problem  which  can  be  modeled  continuously,  the  program 
’SEEK. FOR’  offers  an  easily  used,  interactive  tool  for  the  determination  of  a 
practical  optimum  solution.  The  program  is  small  with  an  object  code  size  of 
46  k  bytes  and  runs  quickly  on  a  486  personal  computer  running  DOS. 
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SUMMARY  OF  RESULTS 


In  this  report,  the  use  of  a  computer  program,  ’SEEK. FOR’  is 
demonstrated  with  four  optimization  examples.  The  program  is  written  to  work 
with  two  user-provided  subroutines  and  an  input  data  file.  Its  purpose  is  to 
perform  a  gradient  search  optimization  of  a  user’s  problem.  The  method  of 
optimization  uses  a  modified  feasible  directions  gradient  in  addition  to 
simpler  gradients  in  the  objective  function  and  the  violated  constraints.  The 
program  is  written  in  ANSI  standard  Fortran  77,  is  about  1,200  lines  long  and 
has  an  object  size  of  about  46  K  bytes  for  the  optimizing  code  for  use  on  a 
personal  computer  running  DOS. 

To  illustrate  the  use  of  the  program,  this  report  documents  four 
optimization  examples:  a  bushing  design,  a  helical  coil  spring  design,  a  gear 
mesh  design  and  a  two-parameter  Weibull  life-reliability  curve  fit.  In  each 
example,  the  theory  of  the  problem  is  described,  the  organization  of  the 
problem  into  an  optimization  framework  is  given,  the  programming  of  the 
modeling  subroutines  is  explained  and  a  numerical  example  is  shown  with  the 
input  and  output  data  files  explained. 

In  the  bushing  design  problem,  two  independent  design  parameters  are 
found  which  minimize  a  single  objective  function  subject  to  three  inequality 
constraints.  In  the  spring  design  problem,  four  design  parameters  are  reduced 
to  two  independent  design  parameters  and  three  separate  objective  functions 
are  minimized  with  three  separate  runs  of  the  same  program.  The  designs 
satisfy  seven  inequality  constraints.  For  the  gear  design  problem,  two  sets 
of  three  independent  design  parameters  are  found  in  two  separate  runs  to 
satisfy  two  opposing  objective  functions.  In  this  model,  fourteen  inequality 
constraints  are  satisfied.  For  the  curve  fitting  example,  a  table  of  x  and  y 
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data  points  is  read  into  the  program  and  two  parameters  are  found  which  fit  a 
curve  to  the  data  with  minimized  dispersion.  One  innocuous  inequality 
constraint  is  used  to  allow  the  program  to  run,  since  the  program  requires  the 
subroutine  to  be  in  place.  As  written,  the  program  accepts  problems  with  up 
to  forty  input  data  constants  and  forty  inequality  constraints.  The  limit  on 
the  number  of  independent  design  parameters  is  fifteen  as  is  the  limit  on  the 
number  of  objective  function  terms. 

This  report  describes  the  use  of  the  optimizing  program,  gives  four 
examples  of  its  use  and  discusses  the  program  and  method  themselves.  The 
program  ’SEEK. FOR’  is  an  adaptable,  analytical  tool  for  finding  optimal 
solutions  to  technical  problems. 
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APPENDIX  A 


SEEK.DOC  FILE 


PROGRAM  SEEK. FOR 

MODIFIED  GRADIENT  OPTIMIZATION  PROGRAM 

WITH  INEQUALITY  CONSTRAINT  AVOIDANCE 

LINK  THIS  PROGRAM  WITH  TWO  PROBLEM  SPECIFIC  ANALYSIS 
SUBROUTINES  AND  RUN  WITH  A  PROBLEM  SPECIFIC  DATA  FILE. 

"NAME". IN  DATA  INPUT  FILE  REQUIRED 

FORMAT: 


LINE  1  :  NTITLE  (50  CHAR.  MAX) 

FOR  THE  CONSTANTS 

LINE  2  :  NCO 

IN  SETS  OF  THREE  LINES  EACH 
LINE  3  -  (lA)  :  CONST(I) 

LINE  4  -  (2A)  :  NCON  (30  CHAR.  MAX) 

LINE  5  -  (3A)  :  NUCON  (12  CHAR.  MAX) 

THEN  FOR  THE  INDEPENDENT  VARIABLES 

:  NX 


LINE  3  +  3*NC0 


LINE  (IB) 
LINE  (2B) 
LINE  (3B) 


IN  SETS  OF  THREE  LINES  EACH 

:  XLOW(I)  ,  XHIGH(I)  , 
:  NVAR  (30  CHAR.  MAX) 

:  NUVAR  (12  CHAR.  MAX) 

THEN  FOR  THE  CONSTRAINTS 


XZ(I) 


LINE  4  +  3*NC0  +  3*NX 

IN  SETS  OF  THREE  LINES 

LINE  (1C) 

LINE  (2C) 

LINE  (3C) 

THEN  FOR  THE  OBJECTIVE 
LINE  5  +  3*NC0  +  3*NX  +  3*NCS 
LINE  6  +  3*NC0  +  3*NX  +  3*NCS 

IN  SETS  OF  THREE  LINES 


LINE  (ID) 
LINE  (2D) 
LINE  (3D) 


:  NCS 
EACH 

:  LIMIT(I)  ,  CSTR(I) 

:  NCSTR  (30  CHAR.  MAX) 

:  NUCSTR  (12  CHAR.  MAX) 
FUNCTION 
:  DIR 
:  NOB 
EACH 

:  WGHTF(I) 

:  NOBJ  (30  CHAR.  MAX) 

:  NUOBJ  (12  CHAR.  MAX) 


TO  HAVE  ACCESS  TO  A  TABLE  OF  X  AND  Y  DATA  VALUES, 

LET  THE  FIRST  CONSTANT  BE  THE  NUMBER  OF  DATA  POINTS 
AND  INCLUDE  THE  WORDS  "DATA  POINTS"  IN  ITS  DESCRIPTION. 

A  DATA  FILE  WITH  X,Y  DATA  PAIRS  AND  THE  NAME  "NAME". DAT 
SHOULD  ALSO  EXIST  IN  THE  SAME  DIRECTORY. 
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UP  TO  200  DATA  PAIRS  WILL  BE  PLACED  IN  THE  COMMON  BLOCK 
COMMON  /CURVE/XDP(200),YDP(200). 

TWO  ANALYSIS  SUBROUTINES  ARE  REQUIRED 

BOUNDS(CONST,NCO,X,NX,VCSTR,NCS)  ;  AND 
VALUES(CONST,NCO,X, NX, OBJECT, NOB) 

BOUNDS  EVALUATES  THE  INEQUALITY  CONSTRAINTS 

VALUES  DETERMINES  THE  OBJECTIVE  FUNCTION’S 
COMPONENT  PROPERTIES 

DEFINITIONS: 

CONST  -  fixed  design  problem  constant 

CSTR  -  CONSTRAINT  LIMIT  VALUE  INCLUDING  DECIMAL  POINT 
DIR  -  OPTIMIZATION  DIRECTION  (  MIN  ,  MAX  ) 

LIMIT  -  CONSTRAINT  LIMIT  BOUND  TYPE  (  UPPER  ,  LOWER  ) 
NCO  -  NUMBER  OF  PROBLEM  CONSTANTS 

NCS  -  NUMBER  OF  INEQUALITY  CONSTRAINTS 

NOB  -  NUMBER  OF  PROPERTIES  IN  MERIT  FUNCTION 

NX  -  NUMBER  OF  INDEPENDENT  VARIABLES 

NCON  -  CONSTANT  NAME 
NUCON  -  CONSTANT  DIMENSION  UNITS 
NVAR  -  VARIABLE  NAME 
NUVAR  -  VARIABLE  DIMENSION  UNITS 
NCSTR  -  CONSTRAINT  NAME 
NUCSTR  -  CONSTRAINT  DIMENSION  UNITS 
NOBJ  -  MERIT  FUNCTION  COMPONENT  NAME 
NUOBJ  -  MERIT  FUNCTION  COMPONENT  DIMENSION  UNITS 
NTITLE  -  DESIGN  PROBLEM  TITLE 
OBJECT  -  MERIT  FUNCTION  COMPONENT  VALUES 
VCSTR  -  CONSTRAINT  FUNCTION  VALUES 
WGHTF  -  WEIGHTING  COEFFICIENT  FOR  COMPONENT  IN 
LINEAR  MERIT  FUNCTION  SUM 
X  -  INDEPENDENT  DESIGN  VARIABLE 

XHIGH  -  HIGH  VARIABLE  VALUE 
XLOW  -  LOW  VARIABLE  VALUE 
XZ  -  INITIAL  VARIABLE  VALUE 

AT  PRESENT,  THE  ARRAY  SIZE  LIMITS  ARE: 

MAX  NCO  =  40, 

MAX  NCS  =  40, 

MAX  NOB  =  15,  AND 
MAX  NX  =15. 
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C  AVOID  THE  USE  OF  SUBROUTINE  AND  COMMON  BLOCK  NAMES 

C  WHICH  ARE  ALREADY  USED  BY  SEEK. 

C 

C  THE  SUBROUTINE  NAMES  USED  BY  SEEK  ARE: 


C  BACK . FSUBRT 

C  BOUNCE . FSUBRT 

C  BOUNDS . extern 

C  CHECK  .  FSUBRT 

C  GRADNT . FSUBRT 

C  MERIT  .  FSUBRT 

C  RESIZE . FSUBRT 

C  SCAN . FSUBRT 

C  SCOUT  ....  .  FSUBRT 

C  SIZE.  . FSUBRT 

C  UNIT . FSUBRT 

C  VALUES . extern 

C  WALL . FSUBRT 

C 

C  THE  COMMON  BLOCK  NAMES  USED  BY  S 

C 

C  CURVE  .  common 

C  PAR . common 

C  VAR . common 

C  UNITS  .  common 


C 

C  TO  PRINT  OUT  INTERMEDIATE  RESULTS  FROM  BOUNDS,  VALUES 

C  OR  ANY  SUBROUTINE  CALLED  BY  BOUNDS  OR  VALUES,  INCLUDE  THE 
C  COMMON  BLOCK  "UNITS"  IN  THAT  SUBROUTINE,  AS  FOLLOWS: 

C 

C  COMMON  /UNITS/NW,NR,NF,ND 
C 

C  WHERE: 

C  NW  =  WRITE  NUMBER  FOR  WRITING  TO  THE  SCREEN 

C  -  WRITE(NW, . . . ) 

C  NR  =  READ  NUMBER  FOR  READING  FROM  THE  KEYBOARD 

C  -  READ(NR,...) 

C  NF  =  WRITE  NUMBER  FOR  WRITING  TO  THE  OUTPUT  FILE 

C  -  WRITE(NF,...)  ,AND 

C  ND  =  READ  NUMBER  FOR  READING  FROM  THE  INPUT  DATA  FILE 

-  READ(ND,...). 
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APPENDIX  B 


SYMBOLS 

Variables 

A  -  AGMA  velocity  factor  constant 

-  independent  design  variable  scaling  constant 
B  -  material  surface  strength  constant  (psi) 

C  -  center  distance  (in) 

C2  -  minimum  gear  tooth  curvature  (in) 

C  -  dynamic  capacity  (lbs) 

Cg  -  gear  dynamic  capacity  (lbs) 

-  tooth  dynamic  capacity  (lbs) 

d^  -  independent  design  variable  scaling  slope 

d. .  -  wire  diameter  (mm) 

w 

D  -  shaft  or  mean  coil  diameter  (mm) 

E  -  elastic  modulus  (psi) 

f  -  face  width  (in) 

F  -  applied  force  (N) 

Vf  -  normalized  feasible  direction  gradient 

G  -  shear  modulus  (MPa) 

h^  -  unloaded  spring  height  (mm) 

Vh  -  normalized  violated  constraint  gradient 

4 

J  -  polar  moment  of  inertia  (mm  )  or  AGMA  bending  stress  factor 

k  -  spring  rate  (N/mm) 

K  -  Wahl  stress  concentration  factor 

w 

i  -  life  (h) 

-  minimum  life  (h) 

0  '  ' 
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i  -  mesh  life  (103 

m  '  n; 

L  -  bearing  length  (mm) 

M  -  objective  function  vector 

\?in  -  normalized  objective  function  gradient 

VM  -  unsealed  objective  function  gradient 

N  -  number  of  active  spring  coils 

d 

Ndes  “  design  factor 

Njj  -  number  of  data  points 

N  -  number  of  inactive  end  coils 

e 

N|r  -  design  factor  in  fatigue 

N  -  number  of  gear  teeth 

g 

Ng  -  static  design  factor 

OD  -  outside  coil  diameter  (mm) 

P  -  contact  pressure  (MPa) 

Pj  -  diametral  pitch  (in”^) 

PV  -  pressure  times  velocity  scoring  factor  (MPa  -  m/s) 
Qy  -  AGMA  surface  quality  factor 

R  -  gear  radius  (in) 

R  -  reliability  ^ 

-  calculated  reliability 
Rq  -  measured  reliability 

S  -  surface  compression  strength  (psi) 

8lC 

Sse  -  shear  endurance  strength  (MPa) 

Ssu  -  shear  ultimate  strength  (MPa) 

S  -  shear  yield  strength  (MPa) 

sy 

Sue  ■  tensile  strength  constant  (MPa) 
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AS  -  step  size 

Tg  -  base  temperature  (°F) 

-  flash  temperature  (°F) 

-  friction  torque  (N-m) 

Tq  -  applied  torque  (lb-in) 

V  -  shear  force  (N) 

3 

V  -T  -  sprinq  coil  cylinder  volume  (mm  ) 

C01  I 

V|^  -  single  violated  constraint 

3 

-  wire  volume  (mm  ) 

V  -  pitch  line  velocity  (ft/sec) 

-  sliding  velocity  (m/s)  or  (ft/sec) 

VV|^  -  single  violated  constraint  gradient 

3 

w  -  weight  density  (kN/m  ) 

W  -  gear  weight  (lbs) 

-  spring  weight  (N) 

X  -  independent  design  parameter  vector 

3 

XD^  -  life  test  data  value  (10  h) 

X|^  -  independent  design  parameter  lower  value 

X^j  -  independent  design  parameter  upper  value 

Xr  -  geometric  temperature  factor 

-  thermal -elastic  temperature  factor 
Jf|-  -  load  sharing  temperature  factor 

V  -  scaled  independent  design  parameter  vector 

AY  -  incremental  step  change  in  scaled  design  parameter  vector 

YD.  -  reliability  data  point  value 

a  -  gear  tooth  involute  angle  (radians) 
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P  -  length  to  diameter  ratio 

6  -  spring  deflection  (mm) 

-  involute  interference  angle  (radians) 

9  -  spring  wire  rotation  angle  (radians) 

//  -  coefficient  of  friction 

u  -  surface  roughness  (RMS) 

V  -  Poisson’s  ratio 

p  -  radius  of  curvature  (in) 

(7jj  -  bending  stress  (psi) 

Oq  -  dispersion  of  data  from  fitted  curve 

Om  -  Hertzian  contact  stress  (psi) 

hi 

Out.  -  Hertzian  contact  stress  at  the  gear  tooth  tip  (psi) 
T  -  shear  stress  (MPa) 

0  -  pressure  angle  (radians) 

u)  -  gear  angular  velocity  (rad/sec) 

Q  -  shaft  speed  (RPM) 


Subscripts 


a  -  alternating 


al  -  pinion  addendum 

a2  -  gear  addendum 

bl  -  pinion  base 

b2  -  gear  base 


d  -  dynamic 

i  -  independent  design  parameter  index 

j  -  optimization  step  index 


141 


k 


-  violated  constraint  index 


m  -  mean 

max  -  maximum 

min  -  minimum 

sol  -  solid  height 

1  -  pinion 

2  -  gear 

10  -  90  percent  reliability 

Superscripts 

a  -  wire  strength  exponent 

b  -  Wei  bull  slope 

p  -  load-life  exponent 
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